Gaussian Process Buckling Constraints in Blade Stiffened Shell Constitutive Subclass

Use Gaussian Processes (statistical machine learning approach) to learn buckling constraints
Created a new subclass TACSGPBladeStiffenedShellConstitutive from the base class TACSBladeStiffenedShellConstitutive, which overrides the buckling constraints and combined buckling-stress failure criterion
Buckling constraints include global + local modes for axial and shear (same as before except more accurate). Also stiffener crippling modes are treated as well.
The TACSBladeStiffenedShellConstitutive used panel length as a DV which was enforced to be equal to the actual panel length from a separate computation of a panel length constraint. In the new class TACSGPBladeStiffenedShellConstitutive, panel width is also required, and both the panel length and width need to be set in at runtime through some additional infrastructure in TACS (not through a constraint => as this doesn't work great with our FUNtoFEM + ESP/CAPS infrastructure).

TODO:

[x] Implement the machine learning buckling constraints (so far just have the closed-form ones implemented).
[x] Update the failure criterion and sensitivity computations using the new buckling constraints.
[x] Add panel width constraint in tacs/constraints
[x] Add internal tests to the TACSGPBladeStiffenedShellConstitutive class for the following subtests:
- [x] Nondimensional parameters
- [x] Critical loads
- [x] GP Model kernel and predictMeanTestData sens tests
[x] Fix so all internal tests pass
- [x] Closed-form case
- [x] Machine learning case
- [x] Add printLevel to each of these tests so you can turn it off
[x] Fix so failure tests pass
- [x] Add ply fraction derivs in DV failure test
- [x] failStrain test
- [x] failDV test
[x] Check the physical correctness of the smeared vs unsmeared stiffness computation when using the energy closed-form
- [x] Double check if the computePanelStiffness() which uses smearedStiffnesses is correct or if we need a new unsmeared routine here.
- [x] Add computeStiffenerCripplingStiffness() routine for panel-like stiffness of the stiffener
[x] Add cython objects for this class to be constructed.
[x] Add full up test of the failure criterion and internal tests from the python level
[ ] Implement the panel length, width constraints as state variables for each TACS component in the FUNtoFEM TACS interface. Probably don't need a new state variable class, just store these values in lists representing each TACS component.
- Compute a chain rule product of panel length, width variables to get the total shape derivatives for these state variables in FUNtoFEM and just update these state variables inside of FUNtoFEM. Don't want to expose these to the optimizer.
[x] Add in the ML data and a way to read in this data using a classmethod in the cython probably.
[x] Question : Should I change to using the smeared stiffener properties for the in plane load computation for global buckling? Answer: the actual in plane loads observed in the panel are what is computed in the local buckling section with B = 0 matrix. Smearing is done for global buckling to use an unstiffened panel closed-form approach. Of course I don't need to do that since I have a stiffened panel closed-form solution directly on the regular properties (and then later I use ML instead).
[ ] Additional tests for the constitutive class
- [x] Add verification test of the buckling constraints (closed-form and ML) implemented in this class. This will make sure that they are giving the correct values (not just that the derivatives are correct).
- [x] Time the constitutive class on the whole datasets and ensure the forward + adjoint eval time is reasonable.
- [x] If the time to compute Ytest on the ML models in each constitutive class don't scale well with 30k+ elements/constitutive objects => then I could register previous forward and adjoint inputs + outputs and return previous computed outputs.
- [ ] Add shape derivative test over the panel length and width constraints
[x] May also need to check whether the axial, shear, stiffener crippling data is duplicated for each element, constitutive object (I would not like to repeat and duplicate the data for each 30k elements if there are ~20k trials for ~100k floats per each trio of GP model objects as that would lead to about 3*10^9 floats of storage). There are ways to prevent an object from being copied, moved, etc. https://www.javatpoint.com/cpp-class-to-prevent-object-copies. [Answer: appears to not be duplicated as I use pointers according to @timryanb]
[ ] Add in sphinx and full documentation of each method
[x] Rename ShearGaussianProcessModel and the other GP classes as TACSShearGaussianProcessModel etc.
[x] Add a train routine feature so that if KS on the GPs are changed, it updates the weights by inverting a matrix (maybe in Cython / Python level)

@sean-engelstad, before you go further down this path, you should review the element routines that already exists. There already exists sub-routines to compute the mass centroid, moments of inertia, etc. Please review them to avoid duplicating code.

Also @sean-engelstad , not sure how you're intending to do the panel length/width computation but you may be able to re-use what I have implemented in the panel length constraint

@A-CGray I decided to replicate what you did with the panel length constraint, with some small modifications for the panel width constraint. I will just add some code into FUNtoFEM to compute chain rule products of panel length, panel width DVs with the coordinate derivatives sensitivities for our workflow.

Updates 5/13/2024:

Finished the prototype of the GaussianProcessModel classes, and the entire TACSGPBladeStiffenedShellConstitutive class in C++ and in Cython. Finished writing the derivatives of the combined failure index using backpropagation this time (which was one thing Ali had wanted to do with the previous subclass).
Had a pretty fatal stack traceback that I debugged with gdb python. Turns out in multiple cases I was calling the superclass destructor which C++ already calls (thus it calls twice).
Fixed a number of issues in setting up the GP model pointers saved in the GP Blade constitutive cython class -> but got that working now. Had to set the default GP model pointers in the constitutive class to NULL and then if axialGP, etc entries are not None I cast them to the .pyx GP model class names and get the pointer to the C++ classes.
Got a working prototype for the tests/constitutive_tests/test_gp_blade_shell_constitutive.py tester which contains a class for testing the ML objects in the workflow and the ML objects not in the workflow (so it switches to closed-form solutions). The derivatives all currently fail and I am working on fixing those now.

Update: internal tests for ML in the loop now pass! Next steps: test the failure criterion, etc. and other external tests on the constitutive object. Below is a finite difference test with h = 1e-8 central difference approximation compared to an analytic derivative implemented by hand.

Real Mode: Internal tests

TACSGPBladeStiffened..testAllTests full results::
        testNondimensionalParameters = 1.8179e-07
        testAxialCriticalLoads = 1.3995e-08
        testShearCriticalLoads = 8.3006e-08
        testStiffenerCripplingLoad = 9.8267e-11
        testAxialGP all tests = 1.6610e-08
        testShearGP all tests = 1.6610e-08
        testCripplingGp all tests = 6.2022e-06
        Overall max rel error = 6.2022e-06

Complex Mode: Internal Tests

TACSGPBladeStiffened..testAllTests full results::
        testNondimensionalParameters = 1.0644e-15
        testAxialCriticalLoads = 1.5943e-16
        testShearCriticalLoads = 1.2744e-13
        testStiffenerCripplingLoad = 0.0000e+00
        testAxialGP all tests = 1.5187e-16
        testShearGP all tests = 1.5832e-16
        testCripplingGp all tests = 2.3355e-15
        Overall max rel error = 1.2744e-13

Complex mode, failStrain test : now passes

Testing the derivative of the failure index w.r.t. strain for object TACSBladeStiffenedShellConstitutive
Max Err: 2.1316e-14 in component 0.
Max REr: 1.0542e-13 in component 4.
          Val[   ]                  Analytic               Approximate                Rel. Error                Abs. Error
d(failure)/de[  0]   -3.8007334173280654e+01   -3.8007334173280633e+01    5.6084654544881153e-16    2.1316282072803006e-14
d(failure)/de[  1]   -3.2745208408550816e+01   -3.2745208408550795e+01    6.5097408472247376e-16    2.1316282072803006e-14
d(failure)/de[  2]    5.9636778818554667e+01    5.9636778818554667e+01    0.0000000000000000e+00    0.0000000000000000e+00
d(failure)/de[  3]    1.2952084768476078e+00    1.2952084768476078e+00    0.0000000000000000e+00    0.0000000000000000e+00
d(failure)/de[  4]   -1.0530991148876019e-03   -1.0530991148874908e-03    1.0542436214503594e-13    1.1102230246251565e-16
d(failure)/de[  5]   -1.0137600641731894e-03   -1.0137600641732471e-03    5.6896653966408030e-14    5.7679555576228836e-17
d(failure)/de[  6]    0.0000000000000000e+00    0.0000000000000000e+00                         -    0.0000000000000000e+00
d(failure)/de[  7]    0.0000000000000000e+00    0.0000000000000000e+00                         -    0.0000000000000000e+00
d(failure)/de[  8]    0.0000000000000000e+00    0.0000000000000000e+00                         -    0.0000000000000000e+00

Current status 5/15/24, complex mode failDV test: some fail Here the DV nums are: [0 - panelLength, 1 - stiffPitch, 2 - panelThick, (3,4,5) - panel ply fracs, 6 - stiffHeight, 7 - stiffThick, (8,9) - stiff ply fracs, 10 - panelWidth]

First without including fails[4] the stiffenerCripplingFailure (so just modes 0-3): all pass

          Val[   ]                  Analytic               Approximate                Rel. Error                Abs. Error
d(failure)/dx[  0]   -2.4917262482262680e-08   -2.4917262482262584e-08    3.8508624743370405e-15    9.5952951056151210e-23
d(failure)/dx[  1]    8.0933741616030596e-06    8.0933741616022278e-06    1.0277374276726334e-13    8.3178635420372293e-19
d(failure)/dx[  2]   -1.9482058118457711e-03   -1.9482058118457898e-03    9.5720263502769913e-15    1.8648277366750676e-17
d(failure)/dx[  3]    6.7302757747794399e-07    6.7302757747499854e-07    4.3764197687246176e-12    2.9454511949584226e-18
d(failure)/dx[  4]    1.1597891871193940e-05    1.1597891871193598e-05    2.9505475175250075e-14    3.4220131069073734e-19
d(failure)/dx[  5]    1.6880611895570803e-06    1.6880611895569016e-06    1.0587528045566319e-13    1.7872395187065737e-19
d(failure)/dx[  6]   -4.5169088937689061e-03   -4.5169088937689009e-03    1.1521530653650334e-15    5.2041704279304213e-18
d(failure)/dx[  7]   -4.5157033072326432e-03   -4.5157033072326406e-03    5.7623033156264800e-16    2.6020852139652106e-18
d(failure)/dx[  8]   -2.5854634340644107e-06   -2.5854634340634307e-06    3.7904891906863540e-13    9.8001711997322549e-19
d(failure)/dx[  9]   -9.5013497223568452e-07   -9.5013497223557398e-07    1.1633904955271950e-13    1.1053779961668619e-19
d(failure)/dx[ 10]    6.4420860130355276e-06    6.4420860130358563e-06    5.1015894986442694e-14    3.2864878353466853e-19

then with including the last failure mode stiffenerCrippling, i.e. fails[4] : some fail namely : 2 - panelThick, 6 - stiffHeight, 7 - stiffThick, (8,9) - stiff ply fracs

Max REr: 1.7883e+00 in component 9.
          Val[   ]                  Analytic               Approximate                Rel. Error                Abs. Error
d(failure)/dx[  0]   -6.2705971332282821e-09   -6.2705971332282821e-09    0.0000000000000000e+00    0.0000000000000000e+00
d(failure)/dx[  1]    1.2217607868101255e-06    1.2217607868103621e-06    1.9360086939627522e-13    2.3653395052076337e-19
d(failure)/dx[  2]   -1.1136900671107634e-03   -1.1135051243354465e-03    1.6609063692208262e-04    1.8494277531687711e-07
d(failure)/dx[  3]    3.0617404435500533e-07    3.0617404435514281e-07    4.4903883195187083e-13    1.3748403525121362e-19
d(failure)/dx[  4]    3.0617404435500533e-07    3.0617404435514281e-07    4.4903883195187083e-13    1.3748403525121362e-19
d(failure)/dx[  5]    3.0617404435500533e-07    3.0617404435514281e-07    4.4903883195187083e-13    1.3748403525121362e-19
d(failure)/dx[  6]   -3.3680506906615644e-03   -3.3679479446752786e-03    3.0506999506418792e-05    1.0274598628585291e-07
d(failure)/dx[  7]   -3.3615820013972115e-03   -3.3608341791145349e-03    2.2251091330950385e-04    7.4782228267657180e-07
d(failure)/dx[  8]   -6.0862821095854106e-07   -2.1827890922525383e-07    1.7883051693760510e+00    3.9034930173328723e-07
d(failure)/dx[  9]   -6.0862821095854127e-07   -2.1827890922525352e-07    1.7883051693760561e+00    3.9034930173328776e-07
d(failure)/dx[ 10]    1.0686834227945663e-06    1.0686834227943772e-06    1.7694679400946116e-13    1.8910010547452255e-19

All the failureDV tests pass now! A few small fixes with the strain transform computation from panelStrain -> stiffenerStrain did the trick. Another problem was a local variable t was supposed to be stiffener thick but was actually panelThick because wasn't redefined locally.

All 10 tests pass now! namely in the test script test_gp_blade_shell_constitutive.py

Complex step test of failure DV test

Max Err: 1.8865e-17 in component 2.
Max REr: 4.3805e-12 in component 3.
          Val[   ]                  Analytic               Approximate                Rel. Error                Abs. Error
d(failure)/dx[  0]   -3.2584854162225586e-08   -3.2584854162225474e-08    3.4524188062079851e-15    1.1249656330721176e-22
d(failure)/dx[  1]    8.0933749272162201e-06    8.0933749272153866e-06    1.0298304818370740e-13    8.3348042009823153e-19
d(failure)/dx[  2]   -1.9474018924833293e-03   -1.9474018924833482e-03    9.6873264188886932e-15    1.8865117801247777e-17
d(failure)/dx[  3]    6.7302764114468927e-07    6.7302764114174107e-07    4.3805096132984838e-12    2.9482040520369991e-18
d(failure)/dx[  4]    1.1597892968325800e-05    1.1597892968325456e-05    2.9651539079076256e-14    3.4389537658524594e-19
d(failure)/dx[  5]    1.6880613492434901e-06    1.6880613492433107e-06    1.0625160433979342e-13    1.7935922658109810e-19
d(failure)/dx[  6]   -4.5232780788233059e-03   -4.5232780788233033e-03    5.7526536476884563e-16    2.6020852139652106e-18
d(failure)/dx[  7]   -4.5084343173967456e-03   -4.5084343173967430e-03    5.7715939299027094e-16    2.6020852139652106e-18
d(failure)/dx[  8]   -3.4248925127256912e-06   -3.4248925127247099e-06    2.8651634051237215e-13    9.8128766939410694e-19
d(failure)/dx[  9]   -6.8563572486820331e-07   -6.8563572486809235e-07    1.6183712730497074e-13    1.1096131609031334e-19
d(failure)/dx[ 10]    6.4420866224409167e-06    6.4420866224412454e-06    5.1015890160465777e-14    3.2864878353466853e-19

Result of all 10 constitutive tests in complex mode

..........
----------------------------------------------------------------------
Ran 10 tests in 0.015s

OK

Fixed a problem with the Newton iteration abs tolerance (too strict and couldn't converge due to floating point limit). Switched to relative tolerance and added max_iter limit for Newton iteration in closed-form solution for shear critical loads.

The class is functional now => next step is to put in the trained ML data and update kernel functions if need be. Then verify that the crit in plane loads with ML surrogate models in TACS make sense (verify the current implementation).

smdogroup / tacs

Gaussian Process Buckling Constraints in Blade Stiffened Shell Constitutive Subclass #311