Possible bug in load transfer scheme

johnjasa commented 4 years ago

Type of issue

What types of issue is it?

Bugfix (non-breaking change which fixes an issue)

Description

The loads transfer scheme used in load_transfer.py is different than that presented in the OpenAeroStruct journal article.

In the transfer loads script, it first computes the moments at the middle of the spar elements but does not transfer the moment to the nodes using the formula of "moment transport" in that the code doesn't consider the moment arm between the node and the middle of the element.

Steps to reproduce issue

Run any aerostructural analysis and examine the load_transfer.py component.

Current behavior

The moment arms are computed from the aerodynamic centers to the structural centers.

Expected behavior

The moment arms may need to be computed from the aerodynamic centers to the structural nodes; this theory should be double-checked as correct before implemented.

If the transfer needs correcting, the compute() and compute_partials() methods within that component would both need updating, as well as all impacted tests, namely any aerostructural test.

kanekosh commented 4 years ago

I have fixed the bug in load_transger.py, though I haven't updated the tests yet so a PR is not ready. I am wondering if I can get quick comments/feedback before updating the tests. The new code can be found in my branch. @johnjasa @shamsheersc19

Verification

To verify the new implementation, I did three checks.

1. Comparison with current implementation

I compared the aerostructural analysis of the new implementation to the current version of mdolab/master. I did not do any aerostructural optimization here.

Here is the comparison of aero and structural responses. The test setup is the same as examples/run_CRM.py except for the mesh sizes.

Mesh size	Version	CL	CD	Tip displacement [m]	Max Von Mises [MPa]
`ny=5, nx=2`	current	0.55486797	0.03976179	0.78911463	55.341826
	new	0.57066144	0.03977402	0.30585531	18.679515
	error	2.85 %	0.0308 %	-61.2 %	-66.24 %
`ny=15, nx=5`	current	0.44880370	0.03592364	0.32427078	27.272997
	new	0.44994858	0.03593000	0.29448746	25.793820
	error	0.255 %	0.0177 %	-9.18 %	-5.42 %
`ny=45, nx=11`	current	0.43865442	0.035831758	0.29876407	28.140761
	new	0.43876558	0.035832641	0.29579435	27.967486
	error	0.0253 %	0.00246 %	-0.994 %	-0.616 %

Since the error for the finest mesh was small enough, I believe the new implementation is correct. Also, I confirmed that the current code has been predicting larger stresses as Sham mentioned, which is good.

2. Partial derivatives

I ran unit tests in transfer, which runs OpenMDAO's check_partials(), and they all passed. Therefore the compute_partials() is consistent with compute(). I also changed the mesh size and it was ok as well.

3. Regression tests

As expected, all aerostructural regression tests have failed. But their values were not far apart: most of them are within a few % errors, except test_aerostruct_wingbox_+weight_analysis and test_aerostruct_wingbox_analysis which had ~20% errors.

Questions

Do you think these are enough to verify the new implementation? If you think of additional tests or comparisons I should perform, I'd appreciate that! I'm a bit hesitant to go ahead and update the tests because there is no safety net there.
In the script you commented the input/output arrays are flattened ones, but they seem not. Furthermore, when you flatten those in compute_partials() you're using ndarray.flatten() with default args, which does its job in a C-style order. Is it okay to fix these comments, or am I missing something?

Thanks!

shamsheersc19 commented 4 years ago

Nice!! Good to see the convergent trends. I’ll have a look at the code and test some optimization cases with the wingbox as well. Thanks!

shamsheersc19 commented 4 years ago

@kanekosh If you have time, you could also do a similar optimization mesh convergence study (for fuelburn and wingbox mass) with this code: https://github.com/mdolab/OpenAeroStruct/blob/master/openaerostruct/docs/wingbox_mpt_opt_example.py (You'll need to make the opt tol smaller than it currently is and then all you need to do is change the numbers of nodes.)

shamsheersc19 commented 4 years ago

I had a look at the compute() code and it looks good to me.
I also tested the new code for 2 wingbox optimization cases (15x5 and 45x7 meshes using the code example above; SLSQP tol=1e-5) and compared them to the older results. Things look good.

Mesh	Version	Fuelburn [kg]	Wingbox mass [kg]
15x5 mesh	Old	97616.1	25199.2
	New	96998.0	24482.2
	% diff	0.6	2.8

45x7 mesh	Old	96967.0	22475.1
	New	96922.0	22434.9
	% diff	0.0	0.2

I am OK with you proceeding with updating the tests and creating a PR. For the tests with large errors, check to see if the errors make sense (probably because of coarse meshes).

johnjasa commented 4 years ago

Shugo, thanks for implementing this, and Sham, thanks for further testing it. The changes look good to me and the output numbers look reasonable, especially given the trends Margarida was seeing when she implemented this herself.

I say go forward with updating the tests and PR. Then we'll have to increment the minor version because this will change the behavior of the code. Get @nwu63 in the mix if needed for that.

kanekosh commented 4 years ago

Thank you guys for the comments!

mdolab / OpenAeroStruct