Bug fix for non-aligned FFD blocks for rotType=0 and generalization of xFraction

[marcomangano]

Purpose

This PR is a combination of a bug-fix and new feature addition that builds on top of a previous, non-merged patch implemented by @Aagaard87 for his WT optimization paper.

Let's look at the bug-fix first. For WT parametrization, the refAxis object of kind rotType=0 is the most suited to handle (child) FFDs which are not aligned with the main system of reference.
In the current code, the local section twist is implemented correctly, but the single-axis scaling parameters are not used in this part of the code.
This PR adds the necessary attributes and includes a pre/post-scaling rotation to ensure that the scaling axes definition is consistent with those of the local FFD section (see a more detailed example in the Testing section below).
This feature is in principle propagated to other rotType as well, but a warning is added as the sensitivities I am obtaining with rotType!=0 are slightly off. This will require more investigation and might be case-dependent, but the current code capabilities are preserved and, anyway, the parametrization cases this bugfix is addressing rely on rotType=0 a-priori.

Following this, another minor but IMHO useful addition is included. Now the user can define yFraction and zFraction alongside xFraction. This means that we can chose the parameterized location of the refAxis node in the FFD section along both the chordwise and "thickness" direction. The motivation for this is briefly illustrated in the Testing section.
Note that the Fraction option in the direction of the FFD-box spanwise axis is not effective, as the reference axis nodes must lie along the local FFD section.
The fix for non-aligned FFDs is also extended to this feature for consistency.

The files affected by this PR are not related to the other major PRs currently open, so there should not be any conflicts.
I tried to add relevant comments where useful, but I am open to add more details based on the reviewers' feedback.

P.S. The explanation of this PR ended up being longer than expected, might probably add a related example in the upcoming MACH tutorial extension.

Type of change

• Bugfix (non-breaking change which fixes an issue)
• New feature (non-breaking change which adds functionality)
• Breaking change (non-backwards-compatible fix or feature)
• Code style update (formatting, renaming)
• Refactoring (no functional changes, no API changes)
• Documentation update
• Maintenance update
• Other (please describe)

Testing

The current code does not propagate scale_x, scale_y, and scale_z attributes to the FFD sections for rotType=0. Newly added test_24_rot0_nonaligned ensures that these scaling attributes are now effective.
Moreover, this test ensures that rotType=0 reference axis can handle (given appropriate input parameters) FFD blocks that are not aligned with the main system of reference, e.g. the blades of a 3-bladed wind turbine rotor. The newly added input parameters rot0ang and rot0axis are used to provide the user control on this.

The test can be used as reference to explain how pyGeo now handles these scaling+rotation operations. We start from an initial "vertical" FFD box which, using the combination of rotType=0, rot0ang=-90, and rot0axis=[1,0,0] for addRefAxis(), is first rotated to have its "spanwise" axis along the y axis. Then, the script scales the 2nd section along the z axis for a "thickness" increase and the 4th section along the x axis for "chord" increase, it adds a +/- 30 deg twist respectively, and finally rotates the deformed FFD back in the initial position. The twist is added to ensure that the operation order is maintained, and the scaling preserves the orthogonality of the FFD in the section plane.
See the initial (blue) and final (red) set of FFD points for test_24_rot0_nonaligned in the picture below:

This is a particular case as the FFD box is already aligned with the main axis and we "flip" the y and z axes, but the same criteria can be applied to a general rotation. In my case specifically, this feature is instrumental to apply the same deformations to all 3 blades of a WT rotor.
Note that for this test I turned off the .ref file training via FD approach because of the relatively higher error with respect to the actual pyGeo outputs.
I suggest to plan a testing infrastructure refactoring to address this and other issues, including an extension of the testing covergage.

Ultimately, the test_23_xyzfraction ensures that the definition of the now-generalized location of the reference axis node in the FFD section is consistent with the user input. This is a new feature, and the current master code would not recognize the additional input arguments.
This feature is again useful in the case of WT analysis and optimization, where the blade chord is not aligned to the x axis.

Checklist

• [NA] I have run flake8 and black to make sure the code adheres to PEP-8 and is consistently formatted
• I have run unit and regression tests which pass locally with my changes
• I have added new tests that prove my fix is effective or that my feature works
• I have added necessary documentation

[marcomangano]

The tests are failing because the Azure yaml file has not been updated as done in #65 following issue #66 about numpy-stl. I would wait for that PRs to be merged into master and consequently into this PR before merging this - just to avoid the same fix on two separate PRs.
Everything else works!

[bbrelje]

I don't really know what's going on with this PR but I'll do my best to review it!

[bbrelje]

Can you comment on why you think the derivative tests are failing for rot cases not equal to 0?

[bbrelje]

I am 99% sure it's fine to do this but only if the rotation is truly implemented as a linear operator as it should be

[marcomangano]

This is a good point. This came directly from the Mads/Sandy fix so I did not 100% thought it through, I will check the rotVbyW function.

[marcomangano]

This is the rotation function ... do you think we should enforce some normalization of the rotation axis?

[bbrelje]

I think there is a possibility that this can lead to a reference axis that lies outside the volume

for example: if you have a typical ffd with the k direction as the ref axis and you rotate it in the yaw axis (so the wing is swept but the inboard edge isn't aligned with flow direction, then you can end up specifying curve ends that lie outside the FFD volume

though to be fair if you really really cambered our existing approach you could end up with a ref axis that lies outside the ffd volume

not sure if this matters or not but figured I'd bring it up

a bigger problem is that I don't think this is backward compatible with the existing default behavior. Right now the y and z frac depend on the xfrac automatically, with the new setup they're fixed at 0.5. I am surprised that this isn't failing reg tests (do we not use the xfrac option anyplace? or possibly our reg tests use rectangular volumes without pretwist or precamber)

[marcomangano]

The example you bring up is interesting... I am actually wondering if our current code could incur in such a problem.
I assume that for more complex FFD blocks you might want to "manually" feed the Ref Axis curve rather than relying on this approach.
However, at least with the current fix, the FFD is rotated back in the "correct" alignment before you define the nodes, similarly to what is done for the scaling fix, so this "yawed" FFD box you bring up should not be an issue. The correct set-up depends on the user, once more.
Anyway, this is not a good excuse for a sloppy fix, and we should take care of any corner case that comes up to mind. Let me sleep on this and get back to you!

About the backwards incompatibility, I am not sure it is actually such a huge problem here.
In the current code the y/z coordinates are given straight up by:

refaxisNodes[place+i,1] = numpy.mean(self.FFD.coef[sectionArr[i+skip,:,:],1])
refaxisNodes[place+i,2] = numpy.mean(self.FFD.coef[sectionArr[i+skip,:,:],2])

So it's a mean over the entire section points rather than the average of min/max. I agree it's a bit brute force, and I am open to improve it.
For a parallelepiped FFD the behaviour is identical, but I see how the behaviour for a cambered airfoil might be a bit different.
Maybe I am missing some major red flag here. Again, let me sleep on this!

[joanibal]

I'm also worried about the changing the default behavior here.
Can we changing it so if only xFraction is given it works the same as before?

[marcomangano]

Yes, I am envisioning something like getting rid of the defaults Fraction=0.5 and instead having an approach similar to before extended to y and z too.
For example:

if xFraction:
       xMin = numpy.min(self.FFD.coef[sectionArr[i+skip,:,:],0])
       xMax = numpy.max(self.FFD.coef[sectionArr[i+skip,:,:],0])
       refaxisNodes[place+i,0] = xFraction*(xMax - xMin) + xMin
else:
       refaxisNodes[place+i,0] = numpy.mean(self.FFD.coef[sectionArr[i+skip,:,:],0])
       
 [same for y and z]

This way the default behavior will be exactly identical

[bbrelje]

@marcomangano I think you are right and I was wrong about the way this piece of code was working

[joanibal]

I'm also worried about the changing the default behavior here.
Can we changing it so if only xFraction is given it works the same as before?

[marcomangano]

Thanks a lot for the thorough review @joanibal ! I getting back at you with a set of commits to address all the comments by Wed

[marcomangano]

I think I have addressed most of your comments!
The main open question is about the sensitivities being off in some cases. I need to dig a bit more into that, in case would you be ok to reformat the fix to only work with rotType=0?

[marcomangano]

Just to give you an idea of the weird derivatives I am getting when switching to different rot types on my WT case:
Reference rotType==0 derivs: I am running successful optimization with these sensitivities

adjoint_derivs:  {'Analysis_L3_V8_TSR780_P0_S10_mx': OrderedDict([
('twist', array([[ -763.5552927 ,   116.12568456, -2155.34420304,  -993.20320918, 4852.29666231,  6245.96972764, 10594.97544573, 13158.47675163, 8713.77256388,   695.6345934 ]])), 
('pitch', array([[120080.31237854]])), 
('chord', array([[-171990.85678235,  146416.65549705,  550734.84415511, 938855.99457252, 1289369.72123263, 1313284.27661121, 1721695.8959731 , 1853893.08719148, 1380207.89218131, 478848.66615274]]))])}
FD_check_last_chord_var:  {'mx': [474707.66964433715], 'fx': []}

Weird rotType==7 derivs:

adjoint_derivs:  {'Analysis_L3_V8_TSR780_P0_S10_mx': OrderedDict([
('twist', array([[  260.6309593 ,  -592.85889657,  -156.80255433,  -939.95130787, -3817.73396433, -4651.38388404, -7083.64933963, -8264.90312159, -5940.96697309, -1511.07434127]])), 
('pitch', array([[3.76870679]])), 
('chord', array([[-112449.52652624,  126351.96413516,  531239.91652338, 970523.83656604, 1332517.43660235, 1360542.8607918 , 1748574.74904779, 1856629.51856273, 1381378.5093449 , 470057.02531971]]))])}
FD_check_last_chord_var:  {'mx': [474666.9928140938], 'fx': []}

Rotations are really off, but also the chord is inconsistent. For the first case, I would not bother too much because my configuration is pretty unconventional and requires rotType=0... but the chord vars should be handled in the same way. Maybe I am overthinking this idk

[joanibal]

Sorry you may have mentioned this somewhere, but just to be clear.

Are the derivatives incorrect in the current version of pyGeo or only in this PR?

[marcomangano]

@joanibal that is a very good question, I actually just mentioned it on slack.
To clarify, the incorrect derivatives are from the current master, see for example below:

working sensitivities with rotType=0 (consistent with Mads' work)

('twist', array([[ -356.04707827,   213.74663014, -2168.80206387, -1228.51426142, 4575.12050082,  6063.3494016 , 10574.06803871, 13155.82446441, 8712.05983958,   587.98126746]])), 
('pitch', array([[119087.22965122]]))])}
FD_pitch_check:  {'mx': [115571.90936543047], 'fx': []}
FD_last_twist_check:  {'mx': [1475.6038246676326], 'fx': []}

weird rotType=7 sensitivities:

('twist', array([[  -69.34707793,  1500.86287043,   526.64606243,  1773.73572847, 7485.5534904 ,  9132.14309081, 14057.49470594, 16475.04198134, 11844.2212678 ,  2720.01802981]])), 
('pitch', array([[196336.67968909]]))])}
FD_pitch_check:   {'mx': [115540.23375306278], 'fx': []}
FD_last_twist_check: {'mx': [1461.197479441762], 'fx': []}

that's why probably I am overthinking it in the current PR

[marcomangano]

Also worth noting, I made another test (test 24 without rotTheta-based twist) to check if the single axis scaling is handled consistently over the different rot types... and it looks like it is, at least in this simple test.
So I assume there is something weird in the way rotTheta is handled in rotType=7?

I guess understanding and fixing that might go beyond this PR, but we should at least take the opportunity to add a warning to prevent people to use these odd derivatives, as for example I did here. Still not satisfied about that because the scaling fix I added actually "could" be used even with other rotTypes (not 100% sure though)

[joanibal]

hmm, well that's not good.

@marcomangano could you open an issue about the derivatives on github?

[marcomangano]

Will do, just need to isolate a bit the problem.
@bbrelje do you think at this point we should keep the warning below?
https://github.com/marcomangano/pygeo/blob/7fda92b3ca58474e1cd22fd6d10d1f98d4e2a7df/pygeo/DVGeometry.py#L1311

This would prevent using the new features included in this PR with problematic rotTypes, but would not backtrack to other potential issues currently on the master branch.
In the best case scenario, these inconsistencies arise for unconventional set-ups like mine. I assume that having derivatives so off would have affected many other optimization cases if it was a wide spread issue.

[bbrelje]

@marcomangano and @joanibal I'm not sure I totally understand the derivative situation and I don't have time right now to really dig into it. Go ahead and open and issue for the broken derivatives.

Are the derivatives broken for all other rot types or just 7? Why would we allow people to use rot0angles with xyz rotations anyway?

Also, on e.g. line 447 and in many other places you use if rot0ang. This will evaluate to FALSE if rot0ang is 0.0 which may or may not be OK depending on the logic there and if it has any importance to initialization. A place where this could be problematic is e.g. if the baseline value is exactly 0.0 at the start of an opt or something. Please double check to make sure this will be ok. generally we do not use variables that can have float values in boolean checks for this reason.

[marcomangano]

Also, on e.g. line 447 and in many other places you use if rot0ang. This will evaluate to FALSE if rot0ang is 0.0 which may or may not be OK depending on the logic there and if it has any importance to initialization. A place where this could be problematic is e.g. if the baseline value is exactly 0.0 at the start of an opt or something. Please double check to make sure this will be ok. generally we do not use variables that can have float values in boolean checks for this reason.

Good point, I will check these if rot0ang and come up with something more elegant/consistent.
I am using this set-up because the default is None and the user would input (maybe I can add more description) a rot0ang only if necessary. If rot0ang=0 the behavior (with a better check) should and would be consistent with rot0ang=None, just adding a set of non-necessary "zero" rotations.
Actually... maybe it makes sense to keep this if as is, so if the user input rotation is zero pyGeo just skips the extra rotations? Anyway, that's a good catch and I was not fully aware this behavior.

Considering optimizations, note that rot0ang is fixed throughout the optimization as it only helps defining the angular offset of the refAxis with respect to the main system of reference, so this is part of the set up similarly to defining the rotation type and the refAxis curve/xFraction.

Are the derivatives broken for all other rot types or just 7? Why would we allow people to use rot0angles with xyz rotations anyway?

About the broken derivatives, I am doing more tests but the biggest mismatches in my cases come from DVs that use rot_theta, which is only used by rotType=[0,7,8]. Let me dig more into that, maybe we can extend the warning to just [7,8].
I agree that if you have a case where rot0ang is needed you should not use rotType=[1,2,3,4,5,6] as rot_theta more easily supersedes rot_x, rot_y, and rot_z, but technically the fixes I made should extend to these other types.

[marcomangano]

FYI, I am not able to investigate more on the sensitivities issue and open a proper issue until Monday - but I am not forgetting about it!

[marcomangano]

I have been debugging a bit the issue with sensitivities for rotType=7 but I could not really understand what is going on. I think the current PR only brings up an underlying issue, but at the same time there is no evidence we can replicate this discrepancy with other cases. I will open a general issue just to warn users to be careful with sensitivities when using that approach.

Regarding this PR, although the original fix by Mads extended to all rotType (not sure if it was ever tested though...), I realized there is no reason to lose our mind over this "generalization", and just keep the working implementation.
For this reason, I restricted the rot0ang and rot0axis to be only effective with rotType=0 (see this commit), in the same way rot_theta only works with rotType=[0,7,8]. If the user needs non-aligned FFD blocks, he should just use rotType=0 right away.
Do we need to add a warning in addrefAxis() to signal the user if one of the provided inputs is not effective? This would work with both my new options and rot_theta when combined with some rotType.

Also, note I changed the check on the input rotation. For completeness, I verified that rot0ang=0 and rot0ang=None provide the same output, as suggested by @nwu63 .

This PR is ready to be merged for me, let me know if you have other concerns! Thx for your thorough review

[ewu63]

Do we need to add a warning in addrefAxis() to signal the user if one of the provided inputs is not effective? This would work with both my new options and rot_theta when combined with some rotType.

IMO if this is documented in the docstrings this is fine.

Also, note I changed the check on the input rotation. For completeness, I verified that rot0ang=0 and rot0ang=None provide the same output, as suggested by @nwu63 .

Good stuff!

[marcomangano]

I opened an issue to keep track of the issue with the sensitivities.
However, after the last commits, there is no way to reproduce the problem with the changes brought by this PR.

[ewu63]

@bbrelje can you take a final look and merge this if all looks good?