BeamDyn finite differencing capability and variational consistent out…

…put in summary file (#155)

* Adding the initial framework for FD. Changing the computations in summary
file to be variationally consistent. Changing the convergence criteria to
be consistent with literature.

* finished implementation of finite difference routines. Adjusted
printing of statements to terminal for a more informed summary of the
simulation.

* Adding the relevant domunation of input file options as well as updated
convergence criteria.
Also fixed passing of error status for the static analysis in BeamDyn

* Update the r-test submodule with new beamdyn input file options

* Bug fix after merge

* Correct a status statement to update “Time”

* Added the finite difference variables to the misc data structure.
Removed excessive write statements for a cleaner output.

* Update r-test submodule

* Re-inserting calls to display simulation time.

docs/source/user/beamdyn/input_files.rst

@@ -45,7 +45,7 @@ or static analysis (``DynamicSolve = false``).
 Gravity Parameters
 ~~~~~~~~~~~~~~~~~~
 
-``Gx`` , ``Gy`` , and ``Gz`` specify the components of gravity vector along :math:`X`, :math:`Y`, and :math:`Z` directions in the global coordinate system, respectively. 
+``Gx`` , ``Gy`` , and ``Gz`` specify the components of gravity vector along :math:`X`, :math:`Y`, and :math:`Z` directions in the global coordinate system, respectively.
 In FAST, this is normally 0, 0, and -9.80665.
 
 Inertial Frame Parameters
@@ -62,7 +62,7 @@ And the following :math:`3 \times 3` direction cosine matrix (``GlbDCM``) relate
    :align: center
 
    Global and blade coordinate systems in BeamDyn.
-   
+
 
 Blade Floating Reference Frame Parameters
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -75,9 +75,9 @@ Based on the driver input file, the floating blade reference fame is assumed to
 
    v_{rt} = \omega_r \times r_t
 
-where :math:`v_{rt}` is the root (origin of the floating blade reference frame) translational velocity vector; :math:`\omega_r` is the constant root (origin of the floating blade reference frame) angular velocity vector; and :math:`r_t` is the global position vector introduced in the previous section at instant :math:`t`, see :numref:`frame`. 
-The floating blade reference frame coincides with the initial floating blade reference frame at the beginning :math:`t=0`.  
-``RootVel(4)``, ``RootVel(5)``, and ``RootVel(6)`` specify the three components of the constant root angular velocity vector about :math:`X`, :math:`Y`, and :math:`Z` axises in global coordinate system, respectively. 
+where :math:`v_{rt}` is the root (origin of the floating blade reference frame) translational velocity vector; :math:`\omega_r` is the constant root (origin of the floating blade reference frame) angular velocity vector; and :math:`r_t` is the global position vector introduced in the previous section at instant :math:`t`, see :numref:`frame`.
+The floating blade reference frame coincides with the initial floating blade reference frame at the beginning :math:`t=0`.
+``RootVel(4)``, ``RootVel(5)``, and ``RootVel(6)`` specify the three components of the constant root angular velocity vector about :math:`X`, :math:`Y`, and :math:`Z` axises in global coordinate system, respectively.
 ``RootVel(1)``, ``RootVel(2)``, and ``RootVel(3)``, which are the three components of the root translational velocity vector along :math:`X`, :math:`Y`, and :math:`Z` directions in global coordinate system, respectively, are calculated based on Eq. :eq:`rootvelocity`.
 
 BeamDyn can handle more complicated root motions by changing, for example, the ``BD_InputSolve`` subroutine in the ``Driver_Beam.f90``
@@ -258,6 +258,12 @@ time-integration in seconds. The keyword “DEFAULT” may be used to
 indicate that the module should employ the time increment prescribed by
 the driver code (FAST/stand-alone driver program).
 
+``load_retries`` specifies the maximum number of load retries allowed. This option
+currently works only for static analysis. For every load retry, the applied load is
+halved to promote convergence of the Newton-Raphson scheme in iteration of smaller
+load steps as opposed to one single large load step which may cause divergence of the
+Newton-Raphson scheme. The keyword “DEFAULT” sets ``load_retries = 20``.
+
 ``NRMax`` specifies the maximum number of iterations per time step in
 the Newton-Raphson scheme. If convergence is not reached within this
 number of iterations, BeamDyn returns an error message and terminates
@@ -270,6 +276,27 @@ iteration. The keyword “DEFAULT” sets
 ``Stop_Tol = 1.0E-05``. Please refer to
 :numref:`convergence-criterion` for more details.
 
+``tngt_stf_fd`` is a boolean that sets the flag to compute the tangent stiffness
+matrix using finite differencing instead of analytical differentiation. The
+finite differencing is performed using a central scheme. The keyword “DEFAULT”
+sets ``tngt_stf_fd = FALSE``.
+
+``tngt_stf_comp`` is a boolean that sets the flag to compare the analytical tangent
+stiffness matrix against the finite differenced tangent stiffness matrix. Information is
+written to the terminal regarding the dof where the maximum difference is observed. If
+``tngt_stf_fd = FALSE`` and ``tngt_stf_comp = TRUE``, the analytical tangent stiffness
+matrix is used to solve the system of equations while the finite difference tangent stiffness
+matrix is used only to perform the comparison of the two matrices. The keyword “DEFAULT”
+sets ``tngt_stf_comp = FALSE``.
+
+``tngt_stf_pert`` sets the perturbation size for finite differencing. The “DEFAULT” value
+based on experience is set to ``1e-06``.
+
+``tngt_stf_difftol`` is the maximum allowable relative difference between the analytical
+and finite differenced tangent stiffness matrices. If for any entry in the matrices,
+the relative difference exceeds this value the simulation will terminate. The “DEFAULT” value
+is currently set to ``1e-01``.
+
 ``RotStates`` is a flag that indicates if BeamDyn's continuous states should be
 oriented in the rotating frame during linearization analysis when coupled to OpenFAST. 
 If multi-blade coordinate (MBC3) analysis is performed, ``RotStates`` must be ``true``.
@@ -311,7 +338,7 @@ the following equality:
 
 .. math::
    :label: keypoint
-           
+
    kp\_total = \sum_{i=1}^{member\_total} n_i - member\_total +1
 
 where :math:`n_i` is the number of key points in the :math:`i^{th}`
@@ -325,16 +352,16 @@ cases for member and key-point definition.
    :width: 60%
    :align: center
 
-   Member and key point definition: one member defined by four key points; 
+   Member and key point definition: one member defined by four key points;
 
 .. _geometry1-case2:
 
 .. figure:: figs/Geometry_Member2.png
    :width: 60%
    :align: center
 
-   Member and key point definition: two members defined by six key points.           
-   
+   Member and key point definition: two members defined by six key points.
+
 The next section defines the key-point information, preceded by two
 header lines. Each key point is defined by three physical coordinates
 (``kp_xr``, ``kp_yr``, ``kp_zr``) in the IEC standard blade
@@ -355,7 +382,7 @@ more details on the blade geometry definition.
 
    BeamDyn Blade Geometry - Top: Side View; Middle: Front View (Looking Downwind); Bottom: Cross Section View (Looking Toward the Tip, from the Root)
 
-   
+
 Mesh Parameter
 ~~~~~~~~~~~~~~
 
@@ -469,7 +496,7 @@ stiffness-proportional damping coefficients, whereby the
 
 .. math::
    :label: damping-force
-           
+
    \mathcal{\underline{F}}^{Damp} = \underline{\underline{\mu}}~\underline{\underline{S}}~\dot{\underline{\epsilon}}
 
 where :math:`\mathcal{\underline{F}}^{Damp}` is the damping force,
@@ -480,8 +507,8 @@ coefficient matrix defined as
 
 .. math::
    :label: damp-matrix
-           
-   \underline{\underline{\mu}} = 
+
+   \underline{\underline{\mu}} =
    \begin{bmatrix}
        \mu_{11} & 0 & 0 & 0 & 0 & 0 \\
        0 & \mu_{22} & 0 & 0 & 0 & 0 \\
@@ -512,7 +539,7 @@ matrix is given as follows:
 
 .. math::
    :label: Stiffness
-           
+
    \begin{bmatrix}
    K_{ShrFlp} & 0 & 0 & 0 & 0 & 0 \\
    0 & K_{ShrEdg} & 0 & 0 & 0 & 0 \\
@@ -557,5 +584,3 @@ We note that for beam structure, the :math:`i_{plr}` is given as
    :label: PolarMOI
 
    i_{plr} = i_{Edg} + i_{Flp}
-
-