add check for valid rotation matrix to beamdyn subroutine BD_CrvExtractCrv

[mjs271]

This commit modifies the subroutine to check whether R is a valid rotation matrix by checking R*R^T = 1 and det(R) = 1. If it is not a valid rotation matrix, the closest rotation matrix (closest orthogonal matrix) is computed by computing the SVD for R = USV^T, and the new matrix to be used is Rr = UV^T.

This commit also adds LAPACK_gesvd to NWTC_LAPACK to handle single/double precision.

The only issue I see at this time, which I discussed with Mike, is the error catching/return in BD_CrvExtractCrv, as the return may only kick control back to the calling subroutine, and neither currently has ErrMsg/ErrStat machinery.

[bjonkman]

Couple of comments/questions based on the git diff:

1. How does this change impact code performance? This routine gets called quite a bit. Having this check in unit tests or debug mode is fine, but I personally would eliminate this extra work/memory in release mode because BeamDyn is already too slow. Besides the extra math, allocating temporary arrays in this routine is going to be time/memory intensive. If this is really the solution, those work arrays should be in the misc vars types, but--more generally--if the problem is indeed that the DCM being used isn't a true DCM, I would rather focus on why that is incorrect.
2. You are correct that the error handling needs some work. Currently if there is an error, the routine calling BD_CrvExtractCrv will not know it, and the return value will be uninitialized.

[michaelasprague]

@bjonkman - thanks for the comments. @mschmidt271 please do a check on performance with/without this feature, and please consider Bonnie's other suggestion around miscvar and improving the error handling.

@bjonkman - we saw this as a robust fix where, if we're given a non-rotation matrix (e.g., when the rest of the code is run in single precision), we find the closest rotation matrix. The check should be quite efficient, but you're right that we need to demonstrate that.

[bjonkman]

@michaelasprague - I was able to fix the problem with the rest of the code being run in single precision by changing the precision of a few more variables in ElastoDyn. It's in my openfast f/Envision-BeamDyn branch (I think I sent you an email about this a few weeks ago). I think Andy found another variable in the BD driver that needed to be double-precision--I will check that you access to this change, too.

[HaymanConsulting]

Should the new subroutine, BD_CheckRotMat(), be added to a module in the NWTC library instead of residing within the BeamDyn module? If so, should an interface be built for various datatypes and should the interface be more generic so that it works for matrices of varying dimensions? Also, we shouldn't be coding ErrStat values to their numeric values but should use the parameters: ErrID_None, ErrID_Info, ErrID_Warn, etc.

[mjs271]

@michaelasprague @bjonkman @ghaymanNREL

After putting some work in on this pull request, I believe I have addressed all the relevant issues.

1. Performance/Using MiscVar: Using Test26 for OpenFAST in double precision, I compared using BD_CheckRotMat to the commit before I started working on this and the difference in run time appeared to be negligible (running ensembles of 10 on my local machine, the average difference was 3 seconds or less on a ~7 minute run, however variance was relatively wide for the run times, leading me to believe this is not a significant difference), so it doesn't appear that keeping the arrays allocated with MiscVar is necessary.

2. After speaking with @michaelasprague, it was decided that, at this time, rather than fixing the rotation matrix and generating a low-level error, an invalid rotation matrix will generate a fatal error (ErrID_Fatal) and stop the simulation. This is documented in the BD_CheckRotMat subroutine, and if, in the future, it should be decided that the fix is desired, then the relevant commented lines can be used to implement the nearest orthogonal matrix fix. As well, it was decided that keeping BD_CheckRotMat in BeamDyn_Subs for the time being is the proper course.

3. Error handling was added to everything in the call tree above BD_CheckRotMat to handle the previously mentioned fatal error, and any error generated by lapack_dgesvd.

4. The main bug in the computation of WM parameters from a rotation matrix appeared to be that the casing structure in BD_CrvExtractCrv was not mutually exclusive between cases, causing it to fall into the first true statement evaluated and possibly generating invalid WM parameters. This is fixed by precomputing the max value and changing the if statements to evaluate based on the computed index. I have verified this using @rafaelmudafort 's feature/unit-test branch and feel fully confident that BD_CrvExtractCrv is now computing the correct WM parameters for an arbitrary rotation matrix.