From 6253e3e99ec91a179de20e344117dee789bb901f Mon Sep 17 00:00:00 2001 From: Derek Slaughter Date: Fri, 16 Aug 2024 14:36:32 +0000 Subject: [PATCH 01/10] BeamDyn, ModMesh, NWTC_Num performance improvements --- modules/beamdyn/src/BeamDyn.f90 | 718 ++++++++++--------- modules/beamdyn/src/BeamDyn_Subs.f90 | 55 +- modules/nwtc-library/src/ModMesh.f90 | 114 +-- modules/nwtc-library/src/ModMesh_Mapping.f90 | 40 +- modules/nwtc-library/src/NWTC_Num.f90 | 8 +- 5 files changed, 499 insertions(+), 436 deletions(-) diff --git a/modules/beamdyn/src/BeamDyn.f90 b/modules/beamdyn/src/BeamDyn.f90 index 2f8dbfd1e5..ad69f87474 100644 --- a/modules/beamdyn/src/BeamDyn.f90 +++ b/modules/beamdyn/src/BeamDyn.f90 @@ -1648,7 +1648,7 @@ subroutine Init_MiscVars( p, u, y, m, ErrStat, ErrMsg ) ! Array for storing the position information for the quadrature points. CALL AllocAry(m%qp%uuu, p%dof_node ,p%nqp,p%elem_total, 'm%qp%uuu displacement at quadrature point',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - CALL AllocAry(m%qp%uup, p%dof_node/2,p%nqp,p%elem_total, 'm%qp%uup displacement prime at quadrature point',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + CALL AllocAry(m%qp%uup, p%dof_node ,p%nqp,p%elem_total, 'm%qp%uup displacement prime at quadrature point',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) CALL AllocAry(m%qp%vvv, p%dof_node ,p%nqp,p%elem_total, 'm%qp%vvv velocity at quadrature point',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) CALL AllocAry(m%qp%vvp, p%dof_node ,p%nqp,p%elem_total, 'm%qp%vvp velocity prime at quadrature point',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) CALL AllocAry(m%qp%aaa, p%dof_node ,p%nqp,p%elem_total, 'm%qp%aaa acceleration at quadrature point',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) @@ -1677,7 +1677,7 @@ subroutine Init_MiscVars( p, u, y, m, ErrStat, ErrMsg ) ! Inertial force terms CALL AllocAry(m%qp%Gi, 6,6, p%nqp,p%elem_total, 'm%qp%Gi gyroscopic at quadrature point',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) CALL AllocAry(m%qp%Ki, 6,6, p%nqp,p%elem_total, 'm%qp%Ki stiffness at quadrature point',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - CALL AllocAry(m%qp%Mi, 6,6, p%nqp,p%elem_total, 'm%qp%Mi mass at quadrature point',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + CALL AllocAry(m%qp%Mi, p%nqp, 6,6, p%elem_total, 'm%qp%Mi mass at quadrature point',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) ! Elastic force terms: \f$ \underline{\underline{\mathcal{O}}} \f$, etc. from equation (19-21) of NREL CP-2C00-60759. CALL AllocAry(m%qp%Oe, 6,6, p%nqp,p%elem_total, 'm%qp%Oe term at quadrature point',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) @@ -2360,48 +2360,44 @@ SUBROUTINE BD_DisplacementQP( nelem, p, x, m ) TYPE(BD_ContinuousStateType), INTENT(IN ) :: x !< Continuous states at t TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< misc/optimization variables + INTEGER(IntKi) :: ErrStat !< index to current element + CHARACTER(ErrMsgLen) :: ErrMsg !< index to current element INTEGER(IntKi) :: idx_qp !< index to the current quadrature point INTEGER(IntKi) :: elem_start !< Node point of first node in current element - INTEGER(IntKi) :: idx_node - CHARACTER(*), PARAMETER :: RoutineName = 'BD_DisplacementQP' - - - DO idx_qp=1,p%nqp - ! Node point before start of this element - elem_start = p%node_elem_idx( nelem,1 ) - - !> ### Calculate the the displacement fields in an element - !! Using equations (27) and (28) \n - !! \f$ \underline{u}\left( \xi \right) = - !! \sum_{i=1}^{p+1} h^i\left( \xi \right) \underline{\hat{u}}^i - !! \f$ \n - !! and \n - !! \f$ \underline{u}^\prime \left( \xi \right) = - !! \sum_{k=1}^{p+1} h^{k\prime} \left( \xi \right) \underline{\hat{u}}^i - !! \f$ - !! - !! | Variable | Value | - !! | :---------: | :------------------------------------------------------------------------- | - !! | \f$ \xi \f$ | Element natural coordinate \f$ \in [-1,1] \f$ | - !! | \f$ k \f$ | Node number of a \f$ p^\text{th} \f$ order Langrangian-interpolant | - !! | \f$ h^i \left( \xi \right ) \f$ | Component of the shape function matrix, \f$ \underline{\underline{N}} \f$ | - !! | \f$ h^{k\prime} \left( \xi \right ) \f$ | \f$ \frac{\mathrm{d}}{\mathrm{d}x_1} h^i \left( \xi \right) \f$ | - !! | \f$ \underline{\hat{u}}^i \f$ | \f$ k^\text{th} \f$ nodal value | + !> ### Calculate the the displacement fields in an element + !! Using equations (27) and (28) \n + !! \f$ \underline{u}\left( \xi \right) = + !! \sum_{i=1}^{p+1} h^i\left( \xi \right) \underline{\hat{u}}^i + !! \f$ \n + !! and \n + !! \f$ \underline{u}^\prime \left( \xi \right) = + !! \sum_{k=1}^{p+1} h^{k\prime} \left( \xi \right) \underline{\hat{u}}^i + !! \f$ + !! + !! | Variable | Value | + !! | :---------: | :------------------------------------------------------------------------- | + !! | \f$ \xi \f$ | Element natural coordinate \f$ \in [-1,1] \f$ | + !! | \f$ k \f$ | Node number of a \f$ p^\text{th} \f$ order Langrangian-interpolant | + !! | \f$ h^i \left( \xi \right ) \f$ | Component of the shape function matrix, \f$ \underline{\underline{N}} \f$ | + !! | \f$ h^{k\prime} \left( \xi \right ) \f$ | \f$ \frac{\mathrm{d}}{\mathrm{d}x_1} h^i \left( \xi \right) \f$ | + !! | \f$ \underline{\hat{u}}^i \f$ | \f$ k^\text{th} \f$ nodal value | + + ! Node point before start of this element + elem_start = p%node_elem_idx(nelem,1) - ! Initialize values for summation - m%qp%uuu(:,idx_qp,nelem) = 0.0_BDKi ! displacement field \f$ \underline{u} \left( \xi \right) \f$ - m%qp%uup(:,idx_qp,nelem) = 0.0_BDKi ! displacement field \f$ \underline{u}^\prime \left( \xi \right) \f$ + ! Use matrix multiplication to interpolate position and position derivative to quadrature points + call LAPACK_DGEMM('N','N', 1.0_BDKi, x%q(1:3,elem_start:elem_start+p%nodes_per_elem-1), p%Shp, 0.0_BDKi, m%qp%uuu(1:3,:,nelem), ErrStat, ErrMsg) + call LAPACK_DGEMM('N','N', 1.0_BDKi, x%q(1:3,elem_start:elem_start+p%nodes_per_elem-1), p%ShpDer, 0.0_BDKi, m%qp%uup(1:3,:,nelem), ErrStat, ErrMsg) - DO idx_node=1,p%nodes_per_elem - m%qp%uuu(1:3,idx_qp,nelem) = m%qp%uuu(1:3,idx_qp,nelem) + p%Shp(idx_node,idx_qp) *x%q(1:3,elem_start - 1 + idx_node) - m%qp%uup(1:3,idx_qp,nelem) = m%qp%uup(1:3,idx_qp,nelem) + p%ShpDer(idx_node,idx_qp)/p%Jacobian(idx_qp,nelem)*x%q(1:3,elem_start - 1 + idx_node) - ENDDO + ! Apply Jacobian to get position derivative with respect to X-axis + do idx_qp = 1, p%nqp + m%qp%uup(1:3,idx_qp,nelem) = m%qp%uup(1:3,idx_qp,nelem) / p%Jacobian(idx_qp,nelem) + end do - !> Calculate \f$ \underline{E}_1 = x_0^\prime + u^\prime \f$ (equation 23). Note E_1 is along the z direction. - m%qp%E1(1:3,idx_qp,nelem) = p%E10(1:3,idx_qp,nelem) + m%qp%uup(1:3,idx_qp,nelem) + !> Calculate \f$ \underline{E}_1 = x_0^\prime + u^\prime \f$ (equation 23). Note E_1 is along the z direction. + m%qp%E1(1:3,:,nelem) = p%E10(1:3,:,nelem) + m%qp%uup(1:3,:,nelem) - ENDDO END SUBROUTINE BD_DisplacementQP @@ -2418,6 +2414,8 @@ SUBROUTINE BD_RotationalInterpQP( nelem, p, x, m ) TYPE(BD_ContinuousStateType), INTENT(IN ) :: x !< Continuous states at t TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< misc/optimization variables + INTEGER(IntKi) :: ErrStat !< index to current element + CHARACTER(ErrMsgLen) :: ErrMsg !< index to current element INTEGER(IntKi) :: idx_qp !< index to the current quadrature point INTEGER(IntKi) :: elem_start !< Node point of first node in current element INTEGER(IntKi) :: idx_node !< index to current GLL point in element @@ -2426,8 +2424,6 @@ SUBROUTINE BD_RotationalInterpQP( nelem, p, x, m ) REAL(BDKi) :: cc(3) REAL(BDKi) :: temp33(3,3) REAL(BDKi) :: DCM_root(3,3) !< DCM for first node - CHARACTER(*), PARAMETER :: RoutineName = 'BD_RotationalInterpQP' - !> ## Calculate the interpolated rotational displacements !! To calculate this, the algorithm given in http://www.nrel.gov/docs/fy14osti/60759.pdf @@ -2482,6 +2478,15 @@ SUBROUTINE BD_RotationalInterpQP( nelem, p, x, m ) ENDDO + ! Use matrix multiplication to interpolate rotation and rotation derivative to quadrature points + ! These rotations do not include the root node rotation at this point (added later in function) + call LAPACK_DGEMM('N','N', 1.0_BDKi, m%Nrrr(:,:,nelem), p%Shp, 0.0_BDKi, m%qp%uuu(4:6,:,nelem), ErrStat, ErrMsg) + call LAPACK_DGEMM('N','N', 1.0_BDKi, m%Nrrr(:,:,nelem), p%ShpDer, 0.0_BDKi, m%qp%uup(4:6,:,nelem), ErrStat, ErrMsg) + + ! Apply Jacobian to get rotation derivative with respect to X-axis + do idx_qp = 1, p%nqp + m%qp%uup(4:6,idx_qp,nelem) = m%qp%uup(4:6,idx_qp,nelem) / p%Jacobian(idx_qp,nelem) + end do ! QP rotational interpolation DO idx_qp=1,p%nqp @@ -2507,16 +2512,9 @@ SUBROUTINE BD_RotationalInterpQP( nelem, p, x, m ) !! | \f$ h^{k\prime} \left( \xi \right ) \f$ | \f$ \frac{\mathrm{d}}{\mathrm{d}x_1} h^i \left( \xi \right) \f$ | !! | \f$ \underline{\hat{r}}^i \f$ | \f$ k^\text{th} \f$ nodal value | - - ! Initialize values for summations - rrr = 0.0_BDKi ! intermediate rotation field for calculation - rrp = 0.0_BDKi - - ! Note: `m%Nrrr` is \f$ \underline{\hat{r}}^i \f$ - DO idx_node=1,p%nodes_per_elem - rrr(1:3) = rrr(1:3) + p%Shp(idx_node,idx_qp) *m%Nrrr(1:3,idx_node,nelem) - rrp(1:3) = rrp(1:3) + p%ShpDer(idx_node,idx_qp)/p%Jacobian(idx_qp,nelem)*m%Nrrr(1:3,idx_node,nelem) - ENDDO + ! Get rotation and rotation derivative at quadrature point (root rotation is not included) + rrr = m%qp%uuu(4:6,idx_qp,nelem) + rrp = m%qp%uup(4:6,idx_qp,nelem) !> **Step 3:** Restore the rigid body rotation at node \f$ \xi \f$ with \n !! \f$ \underline{c}(\xi) = \underline{\hat{c}}^1 \oplus \underline{r}(\xi) \f$ \n @@ -2564,55 +2562,77 @@ SUBROUTINE BD_StifAtDeformedQP( nelem, p, m ) TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< misc/optimization variables INTEGER(IntKi) :: idx_qp !< index counter for quadrature point - INTEGER(IntKi) :: temp_id2 !< Index to last node of previous element - INTEGER(IntKi) :: i,j !< generic counters - REAL(BDKi) :: tempR6(6,6) - REAL(BDKi) :: tempBeta6(6,6) + INTEGER(IntKi) :: idx_Stif0 !< Index to last node of previous element + + ! Initial stiffness matrix index + idx_Stif0 = (nelem-1)*p%nqp + ! Loop through quadrature points + do idx_qp = 1, p%nqp - ! see Bauchau 2011 Flexible Multibody Dynamics p 692-693, section 17.7.2 + ! Initial stiffness matrix index + idx_Stif0 = idx_Stif0 + 1 - ! extract the mass and stiffness matrices for the current element - temp_id2 = (nelem-1)*p%nqp + ! Calculate stiffness and damping matrices for this quadrature point + call Calc_Stif_betaC(m%qp%RR0(:,:,idx_qp,nelem), & + p%Stif0_QP(:,:,idx_Stif0), & + m%qp%Stif(:,:,idx_qp,nelem), & + m%qp%betaC(:,:,idx_qp,nelem)) + end do - DO idx_qp=1,p%nqp +contains + subroutine Calc_Stif_betaC(RR0, Stif0, Stif, betaC) + REAL(BDKi), intent(in) :: RR0(:,:), Stif0(:,:) + REAL(BDKi), intent(inout) :: Stif(:,:), betaC(:,:) + REAL(BDKi) :: tempR6(6,6) + REAL(BDKi) :: tempR6_T(6,6) + REAL(BDKi) :: tempBeta6(6,6) + REAL(BDKi) :: tempBeta_diag(6) + INTEGER(IntKi) :: i, j + + ! see Bauchau 2011 Flexible Multibody Dynamics p 692-693, section 17.7.2 !> RR0 is the rotation tensor at quadrature point \f$ \left(\underline{\underline{R}}\underline{\underline{R}}_0\right) \f$ (3x3) - - ! Setup the temporary matrix for modifying the stiffness matrix. RR0 is changing with time. + + ! Setup the temporary matrix for modifying the stiffness matrix. RR0 is changing with time. tempR6 = 0.0_BDKi - tempBeta6 = 0.0_BDKi - tempR6(1:3,1:3) = m%qp%RR0(:,:,idx_qp,nelem) ! upper left -- translation - tempR6(4:6,4:6) = m%qp%RR0(:,:,idx_qp,nelem) ! lower right -- rotation - !NOTE: Bauchau has the lower right corner multiplied by H - - ! Move damping ratio from material frame to the calculation reference frame - ! This is the following: - ! tempBEta6=matmul(tempR6,matmul(diag(p%beta),transpose(tempR6))) - do j=1,6 - do i=1,6 - ! diagonal of p%beta * TRANSPOSE(tempR6) - tempBeta6(i,j) = p%beta(i)*tempR6(j,i) - enddo - enddo - tempBeta6 = matmul(tempR6,tempBeta6) - - - !> Modify the Mass matrix so it is in the calculation reference frame - !! \f$ \begin{bmatrix} - !! \left(\underline{\underline{R}} \underline{\underline{R}}_0\right) & 0 \\ - !! 0 & \left(\underline{\underline{R}} \underline{\underline{R}}_0\right) - !! \end{bmatrix} - !! \underline{\underline{C}} - !! \begin{bmatrix} - !! \left(\underline{\underline{R}} \underline{\underline{R}}_0\right)^T & 0 \\ - !! 0 & \left(\underline{\underline{R}} \underline{\underline{R}}_0\right)^T - !! \end{bmatrix} \f$ - m%qp%Stif(:,:,idx_qp,nelem) = MATMUL(tempR6,MATMUL(p%Stif0_QP(1:6,1:6,temp_id2+idx_qp),TRANSPOSE(tempR6))) - - ! Now apply the damping - m%qp%betaC(:,:,idx_qp,nelem) = matmul(tempBeta6,m%qp%Stif(:,:,idx_qp,nelem)) - ENDDO + tempR6(1:3,1:3) = RR0 ! upper left -- translation + tempR6(4:6,4:6) = RR0 ! lower right -- rotation + !NOTE: Bauchau has the lower right corner multiplied by H + ! Compute the transpose of tempR6 + tempR6_T = TRANSPOSE(tempR6) + + ! Move damping ratio from material frame to the calculation reference frame + ! This is the following: + ! tempBeta6 = matmul(tempR6, matmul(diag(p%beta), transpose(tempR6))) + + ! Compute tempBeta_diag = beta * tempR6_T (for diagonal elements only) + do j = 1, 6 + tempBeta_diag(j) = p%beta(j) * tempR6_T(j, j) + end do + + ! Compute tempBeta6 using tempBeta_diag + do j = 1, 6 + do i = 1, 6 + tempBeta6(i, j) = tempR6(i, j) * tempBeta_diag(j) + end do + end do + + !> Modify the Mass matrix so it is in the calculation reference frame + !! \f$ \begin{bmatrix} + !! \left(\underline{\underline{R}} \underline{\underline{R}}_0\right) & 0 \\ + !! 0 & \left(\underline{\underline{R}} \underline{\underline{R}}_0\right) + !! \end{bmatrix} + !! \underline{\underline{C}} + !! \begin{bmatrix} + !! \left(\underline{\underline{R}} \underline{\underline{R}}_0\right)^T & 0 \\ + !! 0 & \left(\underline{\underline{R}} \underline{\underline{R}}_0\right)^T + !! \end{bmatrix} \f$ + Stif = matmul(tempR6, matmul(Stif0, tempR6_T)) + + ! Now apply the damping + betaC = matmul(tempBeta6, Stif) + end subroutine END SUBROUTINE BD_StifAtDeformedQP @@ -2630,23 +2650,35 @@ SUBROUTINE BD_QPData_mEta_rho( p, m ) TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< misc/optimization variables INTEGER(IntKi) :: nelem !< index to current element number + INTEGER(IntKi) :: qp_start !< index to start qp indexing for element INTEGER(IntKi) :: idx_qp !< index to the current quadrature point - DO nelem=1,p%elem_total - DO idx_qp=1,p%nqp + do nelem = 1, p%elem_total + qp_start = (nelem-1)*p%nqp + do idx_qp = 1, p%nqp + call Calc_RR0mEta_rho(m%qp%RR0(:,:,idx_qp,nelem), & + p%Mass0_QP(:,:,qp_start+idx_qp), & + m%qp%RR0mEta(:,idx_qp,nelem), & + m%qp%rho(:,:,idx_qp,nelem)) + end do + end do + +contains + subroutine Calc_RR0mEta_rho(RR0, Mass0, RR0mEta, rho) + real(BDKi), intent(in) :: RR0(:,:), Mass0(:,:) + real(BDKi), intent(out) :: RR0mEta(:), rho(:,:) + !> Calculate the new center of mass times mass at the deflected location !! as \f$ \left(\underline{\underline{R}}\underline{\underline{R}}_0\right) m \underline{\eta} \f$ - m%qp%RR0mEta(:,idx_qp,nelem) = MATMUL(m%qp%RR0(:,:,idx_qp,nelem),p%qp%mEta(:,idx_qp,nelem)) + m%qp%RR0mEta(:,idx_qp,nelem) = MATMUL(RR0, p%qp%mEta(:,idx_qp,nelem)) !> Calculate \f$ \rho = \left(\underline{\underline{R}}\underline{\underline{R}}_0\right) !! \underline{\underline{M}}_{2,2} !! \left(\underline{\underline{R}}\underline{\underline{R}}_0\right)^T \f$ where !! \f$ \underline{\underline{M}}_{2,2} \f$ is the inertial terms of the undeflected mass matrix at this quadrature point - m%qp%rho(:,:,idx_qp,nelem) = p%Mass0_QP(4:6,4:6,(nelem-1)*p%nqp+idx_qp) - m%qp%rho(:,:,idx_qp,nelem) = MATMUL(m%qp%RR0(:,:,idx_qp,nelem),MATMUL(m%qp%rho(:,:,idx_qp,nelem),TRANSPOSE(m%qp%RR0(:,:,idx_qp,nelem)))) - ENDDO - ENDDO + rho = MATMUL(RR0, MATMUL(Mass0(4:6,4:6), TRANSPOSE(RR0))) + end subroutine END SUBROUTINE BD_QPData_mEta_rho @@ -2663,99 +2695,117 @@ SUBROUTINE BD_ElasticForce(nelem,p,m,fact) TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< Misc/optimization variables. LOGICAL, INTENT(IN ) :: fact !< Boolean to calculate the Jacobian - REAL(BDKi) :: cet !< for storing the \f$ I_{yy} + I_{zz} \f$ inertia term - REAL(BDKi) :: k1s - REAL(BDKi) :: Wrk33(3,3) - REAL(BDKi) :: tildeE(3,3) - REAL(BDKi) :: C21(3,3) - REAL(BDKi) :: epsi(3,3) - REAL(BDKi) :: mu(3,3) + INTEGER(IntKi) :: idx_qp !< Index to quadrature point currently being calculated if (.not. fact) then - do idx_qp=1,p%nqp - call Calc_Fc_Fd() + call Calc_Fc_Fd(m%qp%RR0(:,:,idx_qp,nelem), & + m%qp%uuu(:,idx_qp,nelem), & + m%qp%E1(:,idx_qp,nelem), & + m%qp%Fc(:,idx_qp,nelem), & + m%qp%Fd(:,idx_qp,nelem)) end do - else - do idx_qp=1,p%nqp - - call Calc_Fc_Fd() - - - !> ###Calculate the \f$ \underline{\underline{\mathcal{O}}} \f$ from equation (19) - !! - !! \f$ \underline{\underline{\mathcal{O}}} = - !! \begin{bmatrix} - !! \underline{\underline{0}} & \underline{\underline{C}}_{11} \tilde{E}_1 - \tilde{F} \\ - !! \underline{\underline{0}} & \underline{\underline{C}}_{21} \tilde{E}_1 - \tilde{M} - !! \end{bmatrix} - !! = \begin{bmatrix} - !! \underline{\underline{0}} & \psi_E - \tilde{F} \\ - !! \underline{\underline{0}} & \mu - \tilde{M} - !! \end{bmatrix} - !! \f$ - Wrk33(:,:) = OuterProduct(m%qp%RR0(1:3,3,idx_qp,nelem), m%qp%RR0(1:3,3,idx_qp,nelem)) ! z-direction in IEC coords - C21(:,:) = m%qp%Stif(4:6,1:3,idx_qp,nelem) + cet*k1s*Wrk33(:,:) - - tildeE = SkewSymMat(m%qp%E1(:,idx_qp,nelem)) - epsi(:,:) = MATMUL(m%qp%Stif(1:3,1:3,idx_qp,nelem),tildeE) ! Stif is RR0 * p%Stif0_QP * RR0^T - mu(:,:) = MATMUL(C21,tildeE) - - m%qp%Oe(:,:,idx_qp,nelem) = 0.0_BDKi - m%qp%Oe(1:3,4:6,idx_qp,nelem) = epsi(1:3,1:3) - SkewSymMat(m%qp%Fc(1:3,idx_qp,nelem)) - m%qp%Oe(4:6,4:6,idx_qp,nelem) = mu(1:3,1:3) - SkewSymMat(m%qp%Fc(4:6,idx_qp,nelem)) - - - !> ###Calculated \f$ \underline{\underline{\mathcal{P}}} \f$ from equation (20) - !! - !! \f$ \underline{\underline{\mathcal{P}}} = - !! \begin{bmatrix} - !! \underline{\underline{0}} & \underline{\underline{0}} \\ - !! \left(\underline{\underline{\bar{C}}}_{11} \tilde{E}_1 \right)^T + \tilde{F} - !! \left(\underline{\underline{\bar{C}}}_{11} \tilde{E}_1 \right)^T - !! \end{bmatrix} - !! = \begin{bmatrix} - !! \underline{\underline{0}} & \underline{\underline{0}} \\ - !! \psi_E^T + \tilde{F} & \mu^T - !! \end{bmatrix} \f$ - m%qp%Pe(:,:,idx_qp,nelem) = 0.0_BDKi - m%qp%Pe(4:6,1:3,idx_qp,nelem) = TRANSPOSE(epsi) + SkewSymMat(m%qp%Fc(1:3,idx_qp,nelem)) - m%qp%Pe(4:6,4:6,idx_qp,nelem) = TRANSPOSE(mu) - - !> ###Calculated \f$ \underline{\underline{\mathcal{Q}}} \f$ from equation (21) - !! - !! \f{eqnarray*}{ - !! \underline{\underline{\mathcal{Q}}} - !! & =& \underline{\underline{\Upsilon}} \underline{\underline{\mathcal{O}}} - !! = \begin{bmatrix} 0 & 0 \\ - !! \tilde{E}_1^T & 0 \end{bmatrix} - !! \underline{\underline{\mathcal{O}}} \\ - !! \begin{bmatrix} 0 & 0 \\ - !! 0 & \underline{\underline{\mathcal{Q}}}_{22} \end{bmatrix} - !! & =& \tilde{E}_1^T \underline{\underline{\mathcal{O}}}_{12} - !! = - \tilde{E}_1 \underline{\underline{\mathcal{O}}}_{12} - !! \f}\n - !! Note: \f$ \tilde{E}_1^T = - \tilde{E}_1 \f$ - m%qp%Qe(:,:,idx_qp,nelem) = 0.0_BDKi - m%qp%Qe(4:6,4:6,idx_qp,nelem) = -MATMUL(tildeE,m%qp%Oe(1:3,4:6,idx_qp,nelem)) + call Calc_Fc_Fd(m%qp%RR0(:,:,idx_qp,nelem), & + m%qp%uuu(:,idx_qp,nelem), & + m%qp%E1(:,idx_qp,nelem), & + m%qp%Fc(:,idx_qp,nelem), & + m%qp%Fd(:,idx_qp,nelem)) + + call Calc_Oe_Pe_Qe(m%qp%RR0(:,:,idx_qp,nelem), & + m%qp%Stif(:,:,idx_qp,nelem), & + m%qp%Fc(:,idx_qp,nelem), & + m%qp%Oe(:,:,idx_qp,nelem), & + m%qp%Pe(:,:,idx_qp,nelem), & + m%qp%Qe(:,:,idx_qp,nelem)) end do - - ENDIF + end if contains - subroutine Calc_Fc_Fd() - REAL(BDKi) :: e1s - REAL(BDKi) :: eee(6) !< intermediate array for calculation Strain and curvature terms of Fc - REAL(BDKi) :: fff(6) !< intermediate array for calculation of the elastic force, Fc - REAL(BDKi) :: R(3,3) !< rotation matrix at quatrature point - REAL(BDKi) :: Rx0p(3) !< \f$ \underline{R} \underline{x}^\prime_0 \f$ - !REAL(BDKi) :: Wrk(3) - + subroutine Calc_Oe_Pe_Qe(RR0, Stif, Fc, Oe, Pe, Qe) + REAL(BDKi), intent(in) :: RR0(:,:), Stif(:,:), Fc(:) + REAL(BDKi), intent(inout) :: Oe(:,:), Pe(:,:), Qe(:,:) + REAL(BDKi) :: Wrk33(3,3) + REAL(BDKi) :: tildeE(3,3) + REAL(BDKi) :: C21(3,3) + REAL(BDKi) :: epsi(3,3) + REAL(BDKi) :: mu(3,3) + REAL(BDKi) :: cet !< for storing the \f$ I_{yy} + I_{zz} \f$ inertia term + REAL(BDKi) :: k1s + + !> ###Calculate the \f$ \underline{\underline{\mathcal{O}}} \f$ from equation (19) + !! + !! \f$ \underline{\underline{\mathcal{O}}} = + !! \begin{bmatrix} + !! \underline{\underline{0}} & \underline{\underline{C}}_{11} \tilde{E}_1 - \tilde{F} \\ + !! \underline{\underline{0}} & \underline{\underline{C}}_{21} \tilde{E}_1 - \tilde{M} + !! \end{bmatrix} + !! = \begin{bmatrix} + !! \underline{\underline{0}} & \psi_E - \tilde{F} \\ + !! \underline{\underline{0}} & \mu - \tilde{M} + !! \end{bmatrix} + !! \f$ + Wrk33 = OuterProduct(RR0(1:3,3), RR0(1:3,3)) ! z-direction in IEC coords + C21 = Stif(4:6,1:3) + cet*k1s*Wrk33(:,:) + + tildeE = SkewSymMat(m%qp%E1(:,idx_qp,nelem)) + epsi = MATMUL(Stif(1:3,1:3),tildeE) ! Stif is RR0 * p%Stif0_QP * RR0^T + mu = MATMUL(C21,tildeE) + + Oe = 0.0_BDKi + Oe(1:3,4:6) = epsi(1:3,1:3) - SkewSymMat(Fc(1:3)) + Oe(4:6,4:6) = mu(1:3,1:3) - SkewSymMat(Fc(4:6)) + + + !> ###Calculated \f$ \underline{\underline{\mathcal{P}}} \f$ from equation (20) + !! + !! \f$ \underline{\underline{\mathcal{P}}} = + !! \begin{bmatrix} + !! \underline{\underline{0}} & \underline{\underline{0}} \\ + !! \left(\underline{\underline{\bar{C}}}_{11} \tilde{E}_1 \right)^T + \tilde{F} + !! \left(\underline{\underline{\bar{C}}}_{11} \tilde{E}_1 \right)^T + !! \end{bmatrix} + !! = \begin{bmatrix} + !! \underline{\underline{0}} & \underline{\underline{0}} \\ + !! \psi_E^T + \tilde{F} & \mu^T + !! \end{bmatrix} \f$ + Pe = 0.0_BDKi + Pe(4:6,1:3) = TRANSPOSE(epsi) + SkewSymMat(Fc(1:3)) + Pe(4:6,4:6) = TRANSPOSE(mu) + + !> ###Calculated \f$ \underline{\underline{\mathcal{Q}}} \f$ from equation (21) + !! + !! \f{eqnarray*}{ + !! \underline{\underline{\mathcal{Q}}} + !! & =& \underline{\underline{\Upsilon}} \underline{\underline{\mathcal{O}}} + !! = \begin{bmatrix} 0 & 0 \\ + !! \tilde{E}_1^T & 0 \end{bmatrix} + !! \underline{\underline{\mathcal{O}}} \\ + !! \begin{bmatrix} 0 & 0 \\ + !! 0 & \underline{\underline{\mathcal{Q}}}_{22} \end{bmatrix} + !! & =& \tilde{E}_1^T \underline{\underline{\mathcal{O}}}_{12} + !! = - \tilde{E}_1 \underline{\underline{\mathcal{O}}}_{12} + !! \f}\n + !! Note: \f$ \tilde{E}_1^T = - \tilde{E}_1 \f$ + Qe(:,:) = 0.0_BDKi + Qe(4:6,4:6) = -MATMUL(tildeE,Oe(1:3,4:6)) + end subroutine + + subroutine Calc_Fc_Fd(RR0, uuu, E1, Fc, Fd) + REAL(BDKi), intent(in) :: RR0(:,:), uuu(:), E1(:) + REAL(BDKi), intent(inout) :: Fc(:), Fd(:) + REAL(BDKi) :: e1s + REAL(BDKi) :: eee(6) !< intermediate array for calculation Strain and curvature terms of Fc + REAL(BDKi) :: fff(6) !< intermediate array for calculation of the elastic force, Fc + REAL(BDKi) :: R(3,3) !< rotation matrix at quatrature point + REAL(BDKi) :: Rx0p(3) !< \f$ \underline{R} \underline{x}^\prime_0 \f$ + REAL(BDKi) :: Wrk(3) + REAL(BDKi) :: cet !< for storing the \f$ I_{yy} + I_{zz} \f$ inertia term + REAL(BDKi) :: k1s !> ### Calculate the 1D strain, \f$ \underline{\epsilon} \f$, equation (5) !! \f$ \underline{\epsilon} = \underline{x}^\prime_0 + \underline{u}^\prime - @@ -2768,9 +2818,9 @@ subroutine Calc_Fc_Fd() !! Note: \f$ \underline{\underline{R}}\underline{\underline{R}}_0 \f$ is used to go from the material basis into the inertial basis !! and the transpose for the other direction. ! eee(1:3) = m%qp%E1(1:3,idx_qp,nelem) - m%qp%RR0(1:3,3,idx_qp,nelem) ! Using RR0 z direction in IEC coords - call BD_CrvMatrixR(m%qp%uuu(4:6,idx_qp,nelem), R) ! Get rotation at QP as a matrix + call BD_CrvMatrixR(uuu(4:6), R) ! Get rotation at QP as a matrix Rx0p = matmul(R,p%E10(:,idx_qp,nelem)) ! Calculate rotated initial tangent - eee(1:3) = m%qp%E1(1:3,idx_qp,nelem) - Rx0p ! Use rotated initial tangent in place of RR0*i1 to eliminate likely mismatch between R0*i1 and x0' + eee(1:3) = E1(1:3) - Rx0p ! Use rotated initial tangent in place of RR0*i1 to eliminate likely mismatch between R0*i1 and x0' !> ### Set the 1D sectional curvature, \f$ \underline{\kappa} \f$, equation (5) !! \f$ \underline{\kappa} = \underline{k} + \underline{\underline{R}}\underline{k}_i \f$ @@ -2835,11 +2885,11 @@ subroutine Calc_Fc_Fd() ! Strain into the material basis (eq (39) of Dymore manual) !Wrk(:) = MATMUL(TRANSPOSE(m%qp%RR0(:,:,idx_qp,nelem)),eee(1:3)) !e1s = Wrk(3) !epsilon_{1} in material basis (for major axis of blade, which is z in the IEC formulation) - e1s = dot_product( m%qp%RR0(:,3,idx_qp,nelem), eee(1:3) ) + e1s = dot_product( RR0(:,3), eee(1:3) ) !Wrk(:) = MATMUL(TRANSPOSE(m%qp%RR0(:,:,idx_qp,nelem)),eee(4:6)) !k1s = Wrk(3) !kappa_{1} in material basis (for major axis of blade, which is z in the IEC formulation) - k1s = dot_product( m%qp%RR0(:,3,idx_qp,nelem), eee(4:6) ) + k1s = dot_product( RR0(:,3), eee(4:6) ) !> Add extension twist coupling terms to the \f$ \underline{F}^c_{a} \f$\n @@ -2858,8 +2908,8 @@ subroutine Calc_Fc_Fd() !! \f$ C_{et} = C_{4,4} + C_{5,5} \f$ ! Refer Section 1.4 in "Dymore User's Manual - Formulation and finite element implementation of beam elements". cet= p%Stif0_QP(4,4,(nelem-1)*p%nqp+idx_qp) + p%Stif0_QP(5,5,(nelem-1)*p%nqp+idx_qp) ! Dymore theory (22) - m%qp%Fc(1:3,idx_qp,nelem) = fff(1:3) + 0.5_BDKi*cet*k1s*k1s*m%qp%RR0(1:3,3,idx_qp,nelem) ! Dymore theory (25a). Note z-axis is the length of blade. - m%qp%Fc(4:6,idx_qp,nelem) = fff(4:6) + cet*e1s*k1s*m%qp%RR0(1:3,3,idx_qp,nelem) ! Dymore theory (25b). Note z-axis is the length of blade. + Fc(1:3) = fff(1:3) + 0.5_BDKi*cet*k1s*k1s*RR0(1:3,3) ! Dymore theory (25a). Note z-axis is the length of blade. + Fc(4:6) = fff(4:6) + cet*e1s*k1s*RR0(1:3,3) ! Dymore theory (25b). Note z-axis is the length of blade. !> ###Calculate \f$ \underline{\mathcal{F}}^d \f$, equation (16) !! \f$ \underline{F}^d = @@ -2870,9 +2920,9 @@ subroutine Calc_Fc_Fd() !! = \begin{bmatrix} \underline{0} \\ !! \left(\underline{\mathcal{F}}^c \times \underline{E}_1 \right)^T !! \end{bmatrix} \f$ - m%qp%Fd(1:3,idx_qp,nelem) = 0.0_BDKi + Fd(1:3) = 0.0_BDKi ! ADP uu0 ref: If E1 is referenced against a different curve than Stif0_QP, there will be strange coupling terms here. - m%qp%Fd(4:6,idx_qp,nelem) = cross_product(m%qp%Fc(1:3,idx_qp,nelem), m%qp%E1(:,idx_qp,nelem)) + Fd(4:6) = cross_product(Fc(1:3), E1(:)) end subroutine Calc_Fc_Fd END SUBROUTINE BD_ElasticForce @@ -2893,32 +2943,29 @@ SUBROUTINE BD_QPDataVelocity( p, x, m ) TYPE(BD_ContinuousStateType), INTENT(IN ) :: x !< Continuous states at t TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< Misc/optimization variables + INTEGER(IntKi) :: ErrStat !< index to current element + CHARACTER(ErrMsgLen) :: ErrMsg !< index to current element INTEGER(IntKi) :: nelem !< index to current element INTEGER(IntKi) :: idx_qp !< index to quadrature point - INTEGER(IntKi) :: idx_node !< index to the GLL node INTEGER(IntKi) :: elem_start !< Starting quadrature point of current element - DO nelem=1,p%elem_total - - elem_start = p%node_elem_idx(nelem,1) - - DO idx_qp=1,p%nqp - - !> Calculate the values for the + ! Calculate the velocity term, velocity prime (derivative of velocity with respect to X-axis), and acceleration terms - ! Initialize to zero for summation - m%qp%vvv(:,idx_qp,nelem) = 0.0_BDKi - m%qp%vvp(:,idx_qp,nelem) = 0.0_BDKi + ! Loop through elements + do nelem = 1, p%elem_total - ! Calculate the velocity term, velocity prime (derivative of velocity with respect to X-axis), and acceleration terms - DO idx_node=1,p%nodes_per_elem - m%qp%vvv(:,idx_qp,nelem) = m%qp%vvv(:,idx_qp,nelem) + p%Shp(idx_node,idx_qp) * x%dqdt(:,elem_start-1+idx_node) - m%qp%vvp(:,idx_qp,nelem) = m%qp%vvp(:,idx_qp,nelem) + p%ShpDer(idx_node,idx_qp)/p%Jacobian(idx_qp,nelem) * x%dqdt(:,elem_start-1+idx_node) - ENDDO + ! Get start index of quadrature points for given element + elem_start = p%node_elem_idx(nelem,1) - ENDDO + ! Use matrix multiplication to interpolate velocity and velocity derivative to quadrature points + call LAPACK_DGEMM('N','N', 1.0_BDKi, x%dqdt(:,elem_start:elem_start+p%nodes_per_elem-1), p%Shp, 0.0_BDKi, m%qp%vvv(:,:,nelem), ErrStat, ErrMsg) + call LAPACK_DGEMM('N','N', 1.0_BDKi, x%dqdt(:,elem_start:elem_start+p%nodes_per_elem-1), p%ShpDer, 0.0_BDKi, m%qp%vvp(:,:,nelem), ErrStat, ErrMsg) - ENDDO + ! Apply Jacobian to get velocity derivative with respect to X-axis + do idx_qp = 1, p%nqp + m%qp%vvp(:,idx_qp,nelem) = m%qp%vvp(:,idx_qp,nelem) / p%Jacobian(idx_qp,nelem) + end do + end do END SUBROUTINE BD_QPDataVelocity @@ -2938,30 +2985,22 @@ SUBROUTINE BD_QPDataAcceleration( p, OtherState, m ) TYPE(BD_OtherStateType), INTENT(IN ) :: OtherState !< Other states at t on input; at t+dt on outputs TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< Misc/optimization variables + INTEGER(IntKi) :: ErrStat !< index to current element + CHARACTER(ErrMsgLen) :: ErrMsg !< index to current element INTEGER(IntKi) :: nelem !< index of current element INTEGER(IntKi) :: idx_qp !< index of current quadrature point INTEGER(IntKi) :: idx_node INTEGER(IntKi) :: elem_start - - - ! Initialize to zero for summation - m%qp%aaa = 0.0_BDKi - - ! Calculate the acceleration term at t+dt (OtherState%acc is at t+dt) - - DO nelem=1,p%elem_total + ! Loop through elements + do nelem = 1, p%elem_total elem_start = p%node_elem_idx(nelem,1) - DO idx_qp=1,p%nqp - DO idx_node=1,p%nodes_per_elem - m%qp%aaa(:,idx_qp,nelem) = m%qp%aaa(:,idx_qp,nelem) + p%Shp(idx_node,idx_qp) * OtherState%acc(:,elem_start-1+idx_node) - END DO - END DO + ! Interpolate the acceleration term at t+dt (OtherState%acc is at t+dt) to quadrature points + call LAPACK_DGEMM('N','N', 1.0_BDKi, OtherState%acc(:,elem_start:elem_start+p%nodes_per_elem-1), p%Shp, 0.0_BDKi, m%qp%aaa(:,:,nelem), ErrStat, ErrMsg) - END DO - + end do END SUBROUTINE BD_QPDataAcceleration @@ -3049,23 +3088,18 @@ SUBROUTINE BD_DissipativeForce( nelem, p, m,fact ) TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< Misc/optimization variables LOGICAL, INTENT(IN ) :: fact - REAL(BDKi) :: SS_ome(3,3) - REAL(BDKi) :: ffd(6) - REAL(BDKi) :: D11(3,3) - REAL(BDKi) :: D12(3,3) - REAL(BDKi) :: D21(3,3) - REAL(BDKi) :: D22(3,3) - REAL(BDKi) :: b11(3,3) - REAL(BDKi) :: b12(3,3) - REAL(BDKi) :: alpha(3,3) - INTEGER(IntKi) :: idx_qp !< index of current quadrature point - - + IF (.NOT. fact) then ! skip all but Fc and Fd terms - DO idx_qp=1,p%nqp - call Calc_FC_FD_ffd() ! this modifies m%qp%Fc and m%qp%Fd + DO idx_qp=1,p%nqp + ! this modifies m%qp%Fc and m%qp%Fd + CALL Calc_FC_FD_ffd(m%qp%E1(:,idx_qp,nelem), & + m%qp%vvv(:,idx_qp,nelem), & + m%qp%vvp(:,idx_qp,nelem), & + m%qp%betaC(:,:,idx_qp,nelem), & + m%qp%Fc(:,idx_qp,nelem), & + m%qp%Fd(:,idx_qp,nelem)) END DO ! bjj: we don't use these values when fact is FALSE, so let's save time and ignore them here, too. @@ -3078,72 +3112,100 @@ SUBROUTINE BD_DissipativeForce( nelem, p, m,fact ) ! m%qp%Yd(:,:,:,nelem) = 0.0_BDKi ELSE -!FIXME: sometime we can condense this with vector arithmetic and removing some variables that aren't needed. DO idx_qp=1,p%nqp - CALL Calc_FC_FD_ffd() ! this sets local variable ffd and modifies m%qp%Fc and m%qp%Fd - - D11 = m%qp%betaC(1:3,1:3,idx_qp,nelem) - D12 = m%qp%betaC(1:3,4:6,idx_qp,nelem) - D21 = m%qp%betaC(4:6,1:3,idx_qp,nelem) - D22 = m%qp%betaC(4:6,4:6,idx_qp,nelem) - - b11(1:3,1:3) = -MATMUL(SkewSymMat(m%qp%E1(:,idx_qp,nelem)),D11) - b12(1:3,1:3) = -MATMUL(SkewSymMat(m%qp%E1(:,idx_qp,nelem)),D12) - - SS_ome = SkewSymMat( m%qp%vvv(4:6,idx_qp,nelem) ) - - ! Compute stiffness matrix Sd - m%qp%Sd(1:3,1:3,idx_qp,nelem) = -MATMUL(D11,SS_ome) - m%qp%Sd(1:3,4:6,idx_qp,nelem) = -MATMUL(D12,SS_ome) - m%qp%Sd(4:6,1:3,idx_qp,nelem) = -MATMUL(D21,SS_ome) - m%qp%Sd(4:6,4:6,idx_qp,nelem) = -MATMUL(D22,SS_ome) - - ! Compute stiffness matrix Pd - m%qp%Pd(:,:,idx_qp,nelem) = 0.0_BDKi - m%qp%Pd(4:6,1:3,idx_qp,nelem) = SkewSymMat(ffd(1:3)) - MATMUL(b11,SS_ome) - m%qp%Pd(4:6,4:6,idx_qp,nelem) = -MATMUL(b12,SS_ome) - - ! Compute stiffness matrix Od - m%qp%Od(:,1:3,idx_qp,nelem) = 0.0_BDKi - alpha = SkewSymMat(m%qp%vvp(1:3,idx_qp,nelem)) - MATMUL(SS_ome,SkewSymMat(m%qp%E1(:,idx_qp,nelem))) - m%qp%Od(1:3,4:6,idx_qp,nelem) = MATMUL(D11,alpha) - SkewSymMat(ffd(1:3)) - m%qp%Od(4:6,4:6,idx_qp,nelem) = MATMUL(D21,alpha) - SkewSymMat(ffd(4:6)) - - ! Compute stiffness matrix Qd - m%qp%Qd(:,:,idx_qp,nelem) = 0.0_BDKi - m%qp%Qd(4:6,4:6,idx_qp,nelem) = -MATMUL(SkewSymMat(m%qp%E1(:,idx_qp,nelem)),m%qp%Od(1:3,4:6,idx_qp,nelem)) - ! Compute gyroscopic matrix Gd - m%qp%Gd(:,1:3,idx_qp,nelem) = 0.0_BDKi - m%qp%Gd(1:3,4:6,idx_qp,nelem) = TRANSPOSE(b11) - m%qp%Gd(4:6,4:6,idx_qp,nelem) = TRANSPOSE(b12) - - ! Compute gyroscopic matrix Xd - m%qp%Xd(:,:,idx_qp,nelem) = 0.0_BDKi - m%qp%Xd(4:6,4:6,idx_qp,nelem) = -MATMUL(SkewSymMat(m%qp%E1(:,idx_qp,nelem)),m%qp%Gd(1:3,4:6,idx_qp,nelem)) - - ! Compute gyroscopic matrix Yd - m%qp%Yd(1:3,:,idx_qp,nelem) = 0.0_BDKi - m%qp%Yd(4:6,1:3,idx_qp,nelem) = b11 - m%qp%Yd(4:6,4:6,idx_qp,nelem) = b12 + ! this sets local variable ffd and modifies m%qp%Fc and m%qp%Fd + CALL Calc_FC_FD_ffd(m%qp%E1(:,idx_qp,nelem), & + m%qp%vvv(:,idx_qp,nelem), & + m%qp%vvp(:,idx_qp,nelem), & + m%qp%betaC(:,:,idx_qp,nelem), & + m%qp%Fc(:,idx_qp,nelem), & + m%qp%Fd(:,idx_qp,nelem)) + + call Calc_Sd_Pd_Od_Qd_Gd_Xd_Yd(m%qp%E1(:,idx_qp,nelem), & + m%qp%vvp(:,idx_qp,nelem), & + m%qp%betaC(:,:,idx_qp,nelem), & + m%qp%Sd(:,:,idx_qp,nelem), & + m%qp%Od(:,:,idx_qp,nelem), & + m%qp%Qd(:,:,idx_qp,nelem), & + m%qp%Gd(:,:,idx_qp,nelem), & + m%qp%Xd(:,:,idx_qp,nelem), & + m%qp%Yd(:,:,idx_qp,nelem), & + m%qp%Pd(:,:,idx_qp,nelem)) END DO ENDIF CONTAINS - SUBROUTINE Calc_FC_FD_ffd() - REAL(BDKi) :: eed(6) - + subroutine Calc_Sd_Pd_Od_Qd_Gd_Xd_Yd(E1, vvp, betaC, Sd, Od, Qd, Gd, Xd, Yd, Pd) + REAL(BDKi), intent(in) :: E1(:), vvp(:), betaC(:,:) + REAL(BDKi), intent(inout) :: Sd(:,:), Od(:,:), Qd(:,:), Gd(:,:), Xd(:,:), Yd(:,:), Pd(:,:) + REAL(BDKi) :: D11(3,3), D12(3,3), D21(3,3), D22(3,3) + REAL(BDKi) :: b11(3,3), b12(3,3) + REAL(BDKi) :: alpha(3,3) + REAL(BDKi) :: SS_ome(3,3) + REAL(BDKi) :: ffd(6) + + D11 = betaC(1:3,1:3) + D12 = betaC(1:3,4:6) + D21 = betaC(4:6,1:3) + D22 = betaC(4:6,4:6) + + b11(1:3,1:3) = -MATMUL(SkewSymMat(E1),D11) + b12(1:3,1:3) = -MATMUL(SkewSymMat(E1),D12) + + SS_ome = SkewSymMat( m%qp%vvv(4:6,idx_qp,nelem) ) + + ! Compute stiffness matrix Sd + Sd(1:3,1:3) = -MATMUL(D11,SS_ome) + Sd(1:3,4:6) = -MATMUL(D12,SS_ome) + Sd(4:6,1:3) = -MATMUL(D21,SS_ome) + Sd(4:6,4:6) = -MATMUL(D22,SS_ome) + + ! Compute stiffness matrix Pd + Pd = 0.0_BDKi + Pd(4:6,1:3) = SkewSymMat(ffd(1:3)) - MATMUL(b11,SS_ome) + Pd(4:6,4:6) = -MATMUL(b12,SS_ome) + + ! Compute stiffness matrix Od + alpha = SkewSymMat(vvp(1:3)) - MATMUL(SS_ome,SkewSymMat(E1)) + Od(:,1:3) = 0.0_BDKi + Od(1:3,4:6) = MATMUL(D11,alpha) - SkewSymMat(ffd(1:3)) + Od(4:6,4:6) = MATMUL(D21,alpha) - SkewSymMat(ffd(4:6)) + + ! Compute stiffness matrix Qd + Qd = 0.0_BDKi + Qd(4:6,4:6) = -MATMUL(SkewSymMat(E1),Od(1:3,4:6)) + + ! Compute gyroscopic matrix Gd + Gd(:,1:3) = 0.0_BDKi + Gd(1:3,4:6) = TRANSPOSE(b11) + Gd(4:6,4:6) = TRANSPOSE(b12) + + ! Compute gyroscopic matrix Xd + Xd = 0.0_BDKi + Xd(4:6,4:6) = -MATMUL(SkewSymMat(E1),Gd(1:3,4:6)) + + ! Compute gyroscopic matrix Yd + Yd(1:3,:) = 0.0_BDKi + Yd(4:6,1:3) = b11 + Yd(4:6,4:6) = b12 + end subroutine + + SUBROUTINE Calc_FC_FD_ffd(E1, vvv, vvp, betaC, Fc, Fd) + REAL(BDKi), intent(in) :: E1(:), vvv(:), vvp(:), betaC(:,:) + REAL(BDKi), intent(inout) :: Fc(:), Fd(:) + REAL(BDKi) :: eed(6), ffd(6) + ! Compute strain rates - eed = m%qp%vvp(1:6,idx_qp,nelem) - eed(1:3) = eed(1:3) + cross_product(m%qp%E1(:,idx_qp,nelem),m%qp%vvv(4:6,idx_qp,nelem)) + eed = vvp + eed(1:3) = eed(1:3) + cross_product(E1,vvv(4:6)) ! Compute dissipative force - ffd(1:6) = MATMUL(m%qp%betaC(:,:,idx_qp,nelem),eed) + ffd(1:6) = MATMUL(betaC(:,:),eed) - m%qp%Fc(1:6,idx_qp,nelem) = m%qp%Fc(1:6,idx_qp,nelem) + ffd - m%qp%Fd(4:6,idx_qp,nelem) = m%qp%Fd(4:6,idx_qp,nelem) + cross_product(ffd(1:3),m%qp%E1(:,idx_qp,nelem)) - + Fc(1:6) = Fc(1:6) + ffd + Fd(4:6) = Fd(4:6) + cross_product(ffd(1:3),E1) END SUBROUTINE Calc_FC_FD_ffd END SUBROUTINE BD_DissipativeForce @@ -3260,21 +3322,21 @@ SUBROUTINE BD_InertialMassMatrix( nelem, p, m ) INTEGER(IntKi) :: i INTEGER(IntKi) :: idx_qp !< index of current quadrature point - do idx_qp=1,p%nqp + m%qp%Mi(:,:,:,nelem) = 0.0_BDKi - m%qp%Mi(:,:,idx_qp,nelem) = 0.0_BDKi + do idx_qp=1,p%nqp ! Set diagonal values for mass DO i=1,3 - m%qp%Mi(i,i,idx_qp,nelem) = p%qp%mmm(idx_qp,nelem) + m%qp%Mi(idx_qp,i,i,nelem) = p%qp%mmm(idx_qp,nelem) ENDDO ! set mass-inertia coupling terms - m%qp%Mi(1:3,4:6,idx_qp,nelem) = -SkewSymMat(m%qp%RR0mEta(:,idx_qp,nelem)) - m%qp%Mi(4:6,1:3,idx_qp,nelem) = SkewSymMat(m%qp%RR0mEta(:,idx_qp,nelem)) + m%qp%Mi(idx_qp,1:3,4:6,nelem) = -SkewSymMat(m%qp%RR0mEta(:,idx_qp,nelem)) + m%qp%Mi(idx_qp,4:6,1:3,nelem) = SkewSymMat(m%qp%RR0mEta(:,idx_qp,nelem)) ! Set inertia terms - m%qp%Mi(4:6,4:6,idx_qp,nelem) = m%qp%rho(:,:,idx_qp,nelem) + m%qp%Mi(idx_qp,4:6,4:6,nelem) = m%qp%rho(:,:,idx_qp,nelem) end do @@ -3771,19 +3833,10 @@ SUBROUTINE Integrate_ElementForce(nelem, p, m) INTEGER(IntKi) :: idx_dof1 CHARACTER(*), PARAMETER :: RoutineName = 'Integrate_ElementForce' - DO i=1,p%nodes_per_elem - DO idx_dof1=1,p%dof_node - - m%elf(idx_dof1,i) = 0.0_BDKi - - DO idx_qp = 1,p%nqp ! dot_product( m%qp%Fc (idx_dof1,:,nelem), p%QPtw_ShpDer( :,i)) - m%elf(idx_dof1,i) = m%elf(idx_dof1,i) - m%qp%Fc (idx_dof1,idx_qp,nelem)*p%QPtw_ShpDer(idx_qp,i) - END DO - - DO idx_qp = 1,p%nqp ! dot_product(m%qp%Ftemp(idx_dof1,:,nelem), p%QPtw_Shp_Jac(:,i,nelem) ) - m%elf(idx_dof1,i) = m%elf(idx_dof1,i) - m%qp%Ftemp(idx_dof1,idx_qp,nelem)*p%QPtw_Shp_Jac(idx_qp,i,nelem) - END DO - + DO i = 1, p%nodes_per_elem + DO idx_dof1 = 1, p%dof_node + m%elf(idx_dof1,i) = -(dot_product(m%qp%Fc(idx_dof1,:,nelem), p%QPtw_ShpDer(:,i)) + & + dot_product(m%qp%Ftemp(idx_dof1,:,nelem), p%QPtw_Shp_Jac(:,i,nelem))) ENDDO ENDDO @@ -3796,31 +3849,28 @@ SUBROUTINE Integrate_ElementMass(nelem, p, m) TYPE(BD_ParameterType), INTENT(IN ) :: p !< Parameters TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< Misc/optimization variables - INTEGER(IntKi) :: idx_qp - INTEGER(IntKi) :: i - INTEGER(IntKi) :: j - INTEGER(IntKi) :: idx_dof1, idx_dof2 CHARACTER(*), PARAMETER :: RoutineName = 'Integrate_ElementMass' - - DO j=1,p%nodes_per_elem - DO idx_dof2=1,p%dof_node - - DO i=1,p%nodes_per_elem - DO idx_dof1=1,p%dof_node - - m%elm(idx_dof1,i,idx_dof2,j) = 0.0_BDKi - - DO idx_qp = 1,p%nqp - m%elm(idx_dof1,i,idx_dof2,j) = m%elm(idx_dof1,i,idx_dof2,j) + m%qp%Mi(idx_dof1,idx_dof2,idx_qp,nelem)*p%QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem) - END DO - - END DO - END DO - + INTEGER(IntKi) :: ErrStat + CHARACTER(ErrMsgLen) :: ErrMsg + INTEGER(IntKi) :: j + INTEGER(IntKi) :: idx_dof2 + ! INTEGER(IntKi) :: idx_qp + ! INTEGER(IntKi) :: i + ! INTEGER(IntKi) :: idx_dof1 + + DO j = 1, p%nodes_per_elem + DO idx_dof2 = 1, p%dof_node + ! DO i = 1, p%nodes_per_elem + ! DO idx_dof1 = 1, p%dof_node + ! do idx_qp = 1, p%nqp + ! m%elm(idx_dof1,i,idx_dof2,j) = m%elm(idx_dof1,i,idx_dof2,j) + (m%qp%Mi(idx_qp,idx_dof1,idx_dof2,nelem),p%QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem)) + ! end do + ! END DO + ! END DO + call LAPACK_gemm('T', 'N', 1.0_R8Ki, m%qp%Mi(:,:,idx_dof2,nelem), p%QPtw_Shp_Shp_Jac(:,:,j,nelem), 0.0_R8Ki, m%elm(:,:,idx_dof2,j), ErrStat, ErrMsg) END DO END DO - END SUBROUTINE Integrate_ElementMass @@ -5582,10 +5632,7 @@ SUBROUTINE BD_CalcForceAcc( u, p, OtherState, m, ErrStat, ErrMsg ) ! Add point forces at GLL points to RHS of equation. - DO j=1,p%node_total - m%RHS(1:3,j) = m%RHS(1:3,j) + m%PointLoadLcl(1:3,j) - m%RHS(4:6,j) = m%RHS(4:6,j) + m%PointLoadLcl(4:6,j) - ENDDO + m%RHS = m%RHS + m%PointLoadLcl ! Now set the root reaction force. @@ -6847,6 +6894,7 @@ subroutine BD_UpdateGlobalRef(u, p, x, OtherState, ErrStat, ErrMsg) character(ErrMsgLen) :: ErrMsg2 ! Temporary Error message real(R8Ki) :: GlbWM_old(3), GlbWM_new(3), GlbWM_diff(3) real(R8Ki) :: GlbRot_old(3, 3), GlbRot_new(3, 3), GlbRot_diff(3, 3) + real(R8Ki) :: NodeRot_old(3) real(R8Ki) :: GlbPos_old(3), GlbPos_new(3) real(R8Ki) :: pos(3), rot(3), trans_vel(3), rot_vel(3), uuN0(3) integer(IntKi) :: i, j, temp_id @@ -6887,8 +6935,8 @@ subroutine BD_UpdateGlobalRef(u, p, x, OtherState, ErrStat, ErrMsg) matmul(GlbRot_new, p%uuN0(1:3, j, i))) ! Update the node orientation rotation of the node - call BD_CrvCompose(x%q(4:6, temp_id), GlbWM_diff, x%q(4:6, temp_id), FLAG_R1R2) - + NodeRot_old = x%q(4:6, temp_id) + call BD_CrvCompose(x%q(4:6, temp_id), GlbWM_diff, NodeRot_old, FLAG_R1R2) end do end do diff --git a/modules/beamdyn/src/BeamDyn_Subs.f90 b/modules/beamdyn/src/BeamDyn_Subs.f90 index 519a40e589..c1b9796a07 100644 --- a/modules/beamdyn/src/BeamDyn_Subs.f90 +++ b/modules/beamdyn/src/BeamDyn_Subs.f90 @@ -242,32 +242,31 @@ SUBROUTINE BD_CrvCompose( rr, pp, qq, flag) REAL(BDKi), INTENT( OUT):: rr(3) !< Composed rotation REAL(BDKi), INTENT(IN ):: pp(3) !< Input rotation 1 REAL(BDKi), INTENT(IN ):: qq(3) !< Input rotation 2 - INTEGER ,INTENT(IN ):: flag !< Option flag + INTEGER, INTENT(IN ):: flag !< Option flag - REAL(BDKi) :: pp0 - REAL(BDKi) :: p(3) - REAL(BDKi) :: qq0 - REAL(BDKi) :: q(3) + REAL(BDKi) :: pp0, p(3) + REAL(BDKi) :: qq0, q(3) REAL(BDKi) :: tr1 - REAL(BDKi) :: Delta1 - REAL(BDKi) :: Delta2 + REAL(BDKi) :: Delta1, Delta2 REAL(BDKi) :: dd1 REAL(BDKi) :: dd2 - ! Set the local values pp and qq allowing for the transpose - - IF(flag==FLAG_R1TR2 .OR. flag==FLAG_R1TR2T) THEN ! "transpose" (negative) of first rotation parameter - p = -pp - ELSE - p = pp - ENDIF - - IF(flag==FLAG_R1R2T .OR. flag==FLAG_R1TR2T) THEN ! "transpose" (negative) of second rotation parameter - q = -qq - ELSE - q = qq - ENDIF + ! Set the local values pp (R1) and qq (R2) and apply transpose based on flag value + select case (flag) + case (FLAG_R1R2) + p = pp ! R1 + q = qq ! R2 + case (FLAG_R1R2T) + p = pp ! R1 + q = -qq ! R2^T + case (FLAG_R1TR2) + p = -pp ! R1^T + q = qq ! R2 + case (FLAG_R1TR2T) + p = -pp ! R1^T + q = -qq ! R2^T + end select !> ## Composing the resulting Wiener-Milenkovic parameter !! @@ -289,7 +288,6 @@ SUBROUTINE BD_CrvCompose( rr, pp, qq, flag) !! !! - ! Calculate pp0 and qq0. See Bauchau for the mathematics here (equations 8 to 9 and interviening text) pp0 = 2.0_BDKi - dot_product(p,p) / 8.0_BDKi ! p_0 @@ -297,19 +295,16 @@ SUBROUTINE BD_CrvCompose( rr, pp, qq, flag) Delta1 = (4.0_BDKi - pp0) * (4.0_BDKi - qq0) ! Delta_1 in Bauchau Delta2 = pp0 * qq0 - dot_product(p,q) ! Delta_2 in Bauchau - dd1 = Delta1 + Delta2 ! Denomimator term for \Delta_2 >= 0 - dd2 = Delta1 - Delta2 ! Denomimator term for \Delta_2 < 0 - ! Rescaling to remove singularities at +/- 2 \pi - ! Note: changed this to test on \Delta_2 (instead of dd1 > dd2) for better consistency with documentation. - IF ( Delta2 >= 0.0_BDKi ) THEN - tr1 = 4.0_BDKi / dd1 + ! Rescaling to remove singularities at +/- 2 \pi + ! Note: changed this to test on \Delta_2 (instead of dd1 > dd2) for better consistency with documentation. + IF (Delta2 >= 0.0_BDKi) THEN + tr1 = 4.0_BDKi / (Delta1 + Delta2) ELSE - tr1 = -4.0_BDKi / dd2 + tr1 = -4.0_BDKi / (Delta1 - Delta2) ENDIF - rr = tr1 * (qq0*p + pp0*q + cross_product(p,q)) - + rr = tr1 * (qq0*p + pp0*q + Cross_Product(p,q)) END SUBROUTINE BD_CrvCompose diff --git a/modules/nwtc-library/src/ModMesh.f90 b/modules/nwtc-library/src/ModMesh.f90 index 6661370353..7ede50f626 100644 --- a/modules/nwtc-library/src/ModMesh.f90 +++ b/modules/nwtc-library/src/ModMesh.f90 @@ -2020,33 +2020,73 @@ SUBROUTINE MeshCopy( SrcMesh, DestMesh, CtrlCode, ErrStat , ErrMess & IF (.NOT. SrcMesh%Initialized) RETURN !bjj: maybe we should first CALL MeshDestroy(DestMesh,ErrStat, ErrMess) - IF ( CtrlCode .EQ. MESH_NEWCOPY .OR. CtrlCode .EQ. MESH_SIBLING .OR. CtrlCode .EQ. MESH_COUSIN ) THEN - - IF (CtrlCode .EQ. MESH_NEWCOPY) THEN - IOS_l = SrcMesh%IOS - Force_l = SrcMesh%FieldMask(MASKID_FORCE) - Moment_l = SrcMesh%FieldMask(MASKID_MOMENT) - Orientation_l = SrcMesh%FieldMask(MASKID_ORIENTATION) - TranslationDisp_l = SrcMesh%FieldMask(MASKID_TRANSLATIONDISP) - TranslationVel_l = SrcMesh%FieldMask(MASKID_TRANSLATIONVEL) - RotationVel_l = SrcMesh%FieldMask(MASKID_ROTATIONVEL) - TranslationAcc_l = SrcMesh%FieldMask(MASKID_TRANSLATIONACC) - RotationAcc_l = SrcMesh%FieldMask(MASKID_ROTATIONACC) - nScalars_l = SrcMesh%nScalars - ELSE ! Sibling or cousin - IOS_l = SrcMesh%IOS ; IF ( PRESENT(IOS) ) IOS_l = IOS - Force_l = .FALSE. ; IF ( PRESENT(Force) ) Force_l = Force - Moment_l = .FALSE. ; IF ( PRESENT(Moment) ) Moment_l = Moment - Orientation_l = .FALSE. ; IF ( PRESENT(Orientation) ) Orientation_l = Orientation - TranslationDisp_l = .FALSE. ; IF ( PRESENT(TranslationDisp) ) TranslationDisp_l = TranslationDisp - TranslationVel_l = .FALSE. ; IF ( PRESENT(TranslationVel) ) TranslationVel_l = TranslationVel - RotationVel_l = .FALSE. ; IF ( PRESENT(RotationVel) ) RotationVel_l = RotationVel - TranslationAcc_l = .FALSE. ; IF ( PRESENT(TranslationAcc) ) TranslationAcc_l = TranslationAcc - RotationAcc_l = .FALSE. ; IF ( PRESENT(RotationAcc) ) RotationAcc_l = RotationAcc - nScalars_l = 0 ; IF ( PRESENT(nScalars) ) nScalars_l = nScalars - END IF - - IF ( CtrlCode .EQ. MESH_NEWCOPY .OR. CtrlCode .EQ. MESH_COUSIN ) THEN + select case (CtrlCode) + case (MESH_NEWCOPY) + IOS_l = SrcMesh%IOS + Force_l = SrcMesh%FieldMask(MASKID_FORCE) + Moment_l = SrcMesh%FieldMask(MASKID_MOMENT) + Orientation_l = SrcMesh%FieldMask(MASKID_ORIENTATION) + TranslationDisp_l = SrcMesh%FieldMask(MASKID_TRANSLATIONDISP) + TranslationVel_l = SrcMesh%FieldMask(MASKID_TRANSLATIONVEL) + RotationVel_l = SrcMesh%FieldMask(MASKID_ROTATIONVEL) + TranslationAcc_l = SrcMesh%FieldMask(MASKID_TRANSLATIONACC) + RotationAcc_l = SrcMesh%FieldMask(MASKID_ROTATIONACC) + nScalars_l = SrcMesh%nScalars + case (MESH_SIBLING, MESH_COUSIN) + IF ( PRESENT(IOS) ) then + IOS_l = IOS + else + IOS_l = SrcMesh%IOS + end if + IF ( PRESENT(Force) ) then + Force_l = Force + else + Force_l = .FALSE. + end if + IF ( PRESENT(Moment) ) then + Moment_l = Moment + else + Moment_l = .FALSE. + end if + IF ( PRESENT(Orientation) ) then + Orientation_l = Orientation + else + Orientation_l = .FALSE. + end if + IF ( PRESENT(TranslationDisp) ) then + TranslationDisp_l = TranslationDisp + else + TranslationDisp_l = .FALSE. + end if + IF ( PRESENT(TranslationVel) ) then + TranslationVel_l = TranslationVel + else + TranslationVel_l = .FALSE. + end if + IF ( PRESENT(RotationVel) ) then + RotationVel_l = RotationVel + else + RotationVel_l = .FALSE. + end if + IF ( PRESENT(TranslationAcc) ) then + TranslationAcc_l = TranslationAcc + else + TranslationAcc_l = .FALSE. + end if + IF ( PRESENT(RotationAcc) ) then + RotationAcc_l = RotationAcc + else + RotationAcc_l = .FALSE. + end if + IF ( PRESENT(nScalars) ) then + nScalars_l = nScalars + else + nScalars_l = 0 + end if + end select + + select case (CtrlCode) + case (MESH_NEWCOPY, MESH_COUSIN) CALL MeshCreate( DestMesh, IOS=IOS_l, Nnodes=SrcMesh%Nnodes, ErrStat=ErrStat, ErrMess=ErrMess & ,Force=Force_l & @@ -2123,7 +2163,7 @@ SUBROUTINE MeshCopy( SrcMesh, DestMesh, CtrlCode, ErrStat , ErrMess & DestMesh%RemapFlag = SrcMesh%RemapFlag - ELSE IF ( CtrlCode .EQ. MESH_SIBLING ) THEN + case (MESH_SIBLING) !bjj: we should make sure the mesh has been committed, otherwise the element lists haven't been created, yet (and thus not shared) IF ( ASSOCIATED(SrcMesh%SiblingMesh) ) THEN ErrStat = ErrID_Fatal @@ -2165,17 +2205,7 @@ SUBROUTINE MeshCopy( SrcMesh, DestMesh, CtrlCode, ErrStat , ErrMess & DestMesh%maxelemlist = SrcMesh%maxelemlist DestMesh%nextelem = SrcMesh%nextelem - - ENDIF - - DO i = 1, NELEMKINDS - IF ( ASSOCIATED(SrcMesh%ElemTable) ) THEN - ENDIF - IF ( ASSOCIATED(DestMesh%ElemTable) ) THEN - ENDIF - ENDDO - - ELSE IF ( CtrlCode .EQ. MESH_UPDATECOPY ) THEN + case (MESH_UPDATECOPY) IF ( SrcMesh%nNodes .NE. DestMesh%nNodes ) THEN ErrStat = ErrID_Fatal @@ -2183,7 +2213,7 @@ SUBROUTINE MeshCopy( SrcMesh, DestMesh, CtrlCode, ErrStat , ErrMess & RETURN ENDIF - ELSE IF ( CtrlCode .EQ. MESH_UPDATEREFERENCE ) THEN + case (MESH_UPDATEREFERENCE) IF ( SrcMesh%nNodes .NE. DestMesh%nNodes ) THEN ErrStat = ErrID_Fatal @@ -2195,11 +2225,11 @@ SUBROUTINE MeshCopy( SrcMesh, DestMesh, CtrlCode, ErrStat , ErrMess & DestMesh%RefOrientation = SrcMesh%RefOrientation DestMesh%RemapFlag = SrcMesh%RemapFlag - ELSE + case default ErrStat = ErrID_Fatal ErrMess = 'MeshCopy: Invalid CtrlCode.' RETURN - ENDIF + end select ! These aren't shared between siblings, so they get copied, no matter what the CtrlCode: diff --git a/modules/nwtc-library/src/ModMesh_Mapping.f90 b/modules/nwtc-library/src/ModMesh_Mapping.f90 index df2723a60b..6d9a7325d1 100644 --- a/modules/nwtc-library/src/ModMesh_Mapping.f90 +++ b/modules/nwtc-library/src/ModMesh_Mapping.f90 @@ -3262,38 +3262,28 @@ FUNCTION GetLoadsScaleFactor( Src ) ! LOCAL: INTEGER :: I, j - REAL(ReKi) :: MaxLoad + REAL(ReKi) :: MaxLoad, MaxForce, MaxMoment + IF ( Src%FIELDMASK( MASKID_FORCE ) ) then + MaxForce = maxval(abs(src%Force)) + else + MaxForce = 0.0_ReKi + end if - GetLoadsScaleFactor = 1.0 - MaxLoad = 0.0 - - IF ( Src%FIELDMASK( MASKID_FORCE ) ) THEN - - DO I=1,Src%Nnodes - DO J=1,3 - MaxLoad = MAX(MaxLoad, ABS(Src%Force(j,I) ) ) - END DO - END DO - - END IF - + IF ( Src%FIELDMASK( MASKID_MOMENT ) ) then + MaxMoment = maxval(abs(src%Moment)) + else + MaxMoment = 0.0_ReKi + end if + + MaxLoad = max(MaxForce, MaxMoment) - IF ( Src%FIELDMASK( MASKID_MOMENT ) ) THEN - - DO I=1,Src%Nnodes - DO J=1,3 - MaxLoad = MAX(MaxLoad, ABS(Src%Moment(j,I) ) ) - END DO - END DO - - END IF - IF ( MaxLoad > 10. ) THEN GetLoadsScaleFactor = 10**MIN( NINT(log10(MaxLoad)), 15 ) ! Let's not get carried away and cause overflow; 10E15 is as far as we'll go + else + GetLoadsScaleFactor = 1.0_ReKi END IF - END FUNCTION GetLoadsScaleFactor !---------------------------------------------------------------------------------------------------------------------------------- SUBROUTINE CreateLoadMap_P_to_P( Src, Dest, MeshMap, ErrStat, ErrMsg ) diff --git a/modules/nwtc-library/src/NWTC_Num.f90 b/modules/nwtc-library/src/NWTC_Num.f90 index 2daf088cca..8d7eb9338d 100644 --- a/modules/nwtc-library/src/NWTC_Num.f90 +++ b/modules/nwtc-library/src/NWTC_Num.f90 @@ -485,7 +485,7 @@ END SUBROUTINE ConvertUnitsToEngr !> This function computes the cross product of two 3-element arrays (resulting in a vector): \n !! cross_product = Vector1 \f$\times\f$ Vector2 \n !! Use cross_product (nwtc_num::cross_product) instead of directly calling a specific routine in the generic interface. - FUNCTION Cross_ProductR4(Vector1, Vector2) result(CProd) + PURE FUNCTION Cross_ProductR4(Vector1, Vector2) result(CProd) ! Argument declarations. @@ -505,7 +505,7 @@ FUNCTION Cross_ProductR4(Vector1, Vector2) result(CProd) END FUNCTION Cross_ProductR4 !======================================================================= !> \copydoc nwtc_num::cross_productr4 - FUNCTION Cross_ProductR4R8(Vector1, Vector2) result(CProd) + PURE FUNCTION Cross_ProductR4R8(Vector1, Vector2) result(CProd) ! Argument declarations. @@ -525,7 +525,7 @@ FUNCTION Cross_ProductR4R8(Vector1, Vector2) result(CProd) END FUNCTION Cross_ProductR4R8 !======================================================================= !> \copydoc nwtc_num::cross_productr4 - FUNCTION Cross_ProductR8(Vector1, Vector2) result(CProd) + PURE FUNCTION Cross_ProductR8(Vector1, Vector2) result(CProd) ! Argument declarations. @@ -545,7 +545,7 @@ FUNCTION Cross_ProductR8(Vector1, Vector2) result(CProd) END FUNCTION Cross_ProductR8 !======================================================================= !> \copydoc nwtc_num::cross_productr4 - FUNCTION Cross_ProductR8R4(Vector1, Vector2) result(CProd) + PURE FUNCTION Cross_ProductR8R4(Vector1, Vector2) result(CProd) ! Argument declarations. From c8b984824924d971e5c1d20c5886d62bddbd342c Mon Sep 17 00:00:00 2001 From: Derek Slaughter Date: Mon, 19 Aug 2024 20:18:56 +0000 Subject: [PATCH 02/10] Fix bug in BeamDyn performance commit --- modules/beamdyn/src/BeamDyn.f90 | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/modules/beamdyn/src/BeamDyn.f90 b/modules/beamdyn/src/BeamDyn.f90 index ad69f87474..d37909b0e8 100644 --- a/modules/beamdyn/src/BeamDyn.f90 +++ b/modules/beamdyn/src/BeamDyn.f90 @@ -2695,9 +2695,9 @@ SUBROUTINE BD_ElasticForce(nelem,p,m,fact) TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< Misc/optimization variables. LOGICAL, INTENT(IN ) :: fact !< Boolean to calculate the Jacobian - - - INTEGER(IntKi) :: idx_qp !< Index to quadrature point currently being calculated + REAL(BDKi) :: cet !< for storing the \f$ I_{yy} + I_{zz} \f$ inertia term + REAL(BDKi) :: k1s + INTEGER(IntKi) :: idx_qp !< Index to quadrature point currently being calculated if (.not. fact) then @@ -2734,8 +2734,6 @@ subroutine Calc_Oe_Pe_Qe(RR0, Stif, Fc, Oe, Pe, Qe) REAL(BDKi) :: C21(3,3) REAL(BDKi) :: epsi(3,3) REAL(BDKi) :: mu(3,3) - REAL(BDKi) :: cet !< for storing the \f$ I_{yy} + I_{zz} \f$ inertia term - REAL(BDKi) :: k1s !> ###Calculate the \f$ \underline{\underline{\mathcal{O}}} \f$ from equation (19) !! @@ -2804,8 +2802,6 @@ subroutine Calc_Fc_Fd(RR0, uuu, E1, Fc, Fd) REAL(BDKi) :: R(3,3) !< rotation matrix at quatrature point REAL(BDKi) :: Rx0p(3) !< \f$ \underline{R} \underline{x}^\prime_0 \f$ REAL(BDKi) :: Wrk(3) - REAL(BDKi) :: cet !< for storing the \f$ I_{yy} + I_{zz} \f$ inertia term - REAL(BDKi) :: k1s !> ### Calculate the 1D strain, \f$ \underline{\epsilon} \f$, equation (5) !! \f$ \underline{\epsilon} = \underline{x}^\prime_0 + \underline{u}^\prime - From 1d313c9370587c6fa32088eddea2be64f3a16777 Mon Sep 17 00:00:00 2001 From: Derek Slaughter Date: Tue, 20 Aug 2024 20:54:00 +0000 Subject: [PATCH 03/10] Fix beamdyn_utest incorrect size of m%qp%upp --- modules/beamdyn/tests/test_tools.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/modules/beamdyn/tests/test_tools.F90 b/modules/beamdyn/tests/test_tools.F90 index 1f64ec584e..4850f1770b 100644 --- a/modules/beamdyn/tests/test_tools.F90 +++ b/modules/beamdyn/tests/test_tools.F90 @@ -173,8 +173,8 @@ type(BD_MiscVarType) function simpleMiscVarType(nqp, dof_node, elem_total, nodes call AllocAry(m%qp%RR0mEta, 3, nqp, elem_total, 'qp_RR0mEta', ErrStat, ErrMsg) call AllocAry(m%DistrLoad_QP, 6, nqp, elem_total, 'DistrLoad_QP', ErrStat, ErrMsg) - CALL AllocAry(m%qp%uuu, dof_node ,nqp,elem_total, 'm%qp%uuu displacement at quadrature point',ErrStat,ErrMsg) - CALL AllocAry(m%qp%uup, dof_node/2,nqp,elem_total, 'm%qp%uup displacement prime at quadrature point',ErrStat,ErrMsg) + call AllocAry(m%qp%uuu, dof_node, nqp, elem_total, 'm%qp%uuu displacement at quadrature point', ErrStat, ErrMsg) + call AllocAry(m%qp%uup, dof_node, nqp, elem_total, 'm%qp%uup displacement prime at quadrature point', ErrStat, ErrMsg) ! E1, kappa -- used in force calculations CALL AllocAry(m%qp%E1, dof_node/2,nqp,elem_total, 'm%qp%E1 at quadrature point',ErrStat,ErrMsg) From 48a9e87d1fa3d0040a50c118457bd3ca7ee53380 Mon Sep 17 00:00:00 2001 From: Derek Slaughter Date: Thu, 22 Aug 2024 12:56:48 +0000 Subject: [PATCH 04/10] Fix more bugs in BeamDyn performance improvements --- modules/beamdyn/src/BeamDyn.f90 | 82 +++++++++++++++++++-------------- 1 file changed, 48 insertions(+), 34 deletions(-) diff --git a/modules/beamdyn/src/BeamDyn.f90 b/modules/beamdyn/src/BeamDyn.f90 index d37909b0e8..172fca1bd9 100644 --- a/modules/beamdyn/src/BeamDyn.f90 +++ b/modules/beamdyn/src/BeamDyn.f90 @@ -2606,17 +2606,16 @@ subroutine Calc_Stif_betaC(RR0, Stif0, Stif, betaC) ! This is the following: ! tempBeta6 = matmul(tempR6, matmul(diag(p%beta), transpose(tempR6))) - ! Compute tempBeta_diag = beta * tempR6_T (for diagonal elements only) - do j = 1, 6 - tempBeta_diag(j) = p%beta(j) * tempR6_T(j, j) - end do - - ! Compute tempBeta6 using tempBeta_diag - do j = 1, 6 - do i = 1, 6 - tempBeta6(i, j) = tempR6(i, j) * tempBeta_diag(j) - end do - end do + ! Move damping ratio from material frame to the calculation reference frame + ! This is the following: + ! tempBEta6=matmul(tempR6,matmul(diag(p%beta),transpose(tempR6))) + do j=1,6 + do i=1,6 + ! diagonal of p%beta * TRANSPOSE(tempR6) + tempBeta6(i,j) = p%beta(i)*tempR6(j,i) + enddo + enddo + tempBeta6 = matmul(tempR6,tempBeta6) !> Modify the Mass matrix so it is in the calculation reference frame !! \f$ \begin{bmatrix} @@ -2695,8 +2694,8 @@ SUBROUTINE BD_ElasticForce(nelem,p,m,fact) TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< Misc/optimization variables. LOGICAL, INTENT(IN ) :: fact !< Boolean to calculate the Jacobian - REAL(BDKi) :: cet !< for storing the \f$ I_{yy} + I_{zz} \f$ inertia term - REAL(BDKi) :: k1s + REAL(BDKi) :: cet_t !< for storing the \f$ I_{yy} + I_{zz} \f$ inertia term + REAL(BDKi) :: k1s_t INTEGER(IntKi) :: idx_qp !< Index to quadrature point currently being calculated @@ -2704,21 +2703,33 @@ SUBROUTINE BD_ElasticForce(nelem,p,m,fact) do idx_qp=1,p%nqp call Calc_Fc_Fd(m%qp%RR0(:,:,idx_qp,nelem), & m%qp%uuu(:,idx_qp,nelem), & + p%E10(:,idx_qp,nelem), & m%qp%E1(:,idx_qp,nelem), & + m%qp%kappa(1:3,idx_qp,nelem), & + p%Stif0_QP(:,:,(nelem-1)*p%nqp+idx_qp), & + m%qp%Stif(:,:,idx_qp,nelem), & m%qp%Fc(:,idx_qp,nelem), & - m%qp%Fd(:,idx_qp,nelem)) + m%qp%Fd(:,idx_qp,nelem), & + cet_t, k1s_t) end do else do idx_qp=1,p%nqp call Calc_Fc_Fd(m%qp%RR0(:,:,idx_qp,nelem), & m%qp%uuu(:,idx_qp,nelem), & + p%E10(:,idx_qp,nelem), & m%qp%E1(:,idx_qp,nelem), & + m%qp%kappa(1:3,idx_qp,nelem), & + p%Stif0_QP(:,:,(nelem-1)*p%nqp+idx_qp), & + m%qp%Stif(:,:,idx_qp,nelem), & m%qp%Fc(:,idx_qp,nelem), & - m%qp%Fd(:,idx_qp,nelem)) + m%qp%Fd(:,idx_qp,nelem), & + cet_t, k1s_t) call Calc_Oe_Pe_Qe(m%qp%RR0(:,:,idx_qp,nelem), & m%qp%Stif(:,:,idx_qp,nelem), & + m%qp%E1(:,idx_qp,nelem), & m%qp%Fc(:,idx_qp,nelem), & + cet_t, k1s_t, & m%qp%Oe(:,:,idx_qp,nelem), & m%qp%Pe(:,:,idx_qp,nelem), & m%qp%Qe(:,:,idx_qp,nelem)) @@ -2726,8 +2737,8 @@ SUBROUTINE BD_ElasticForce(nelem,p,m,fact) end if contains - subroutine Calc_Oe_Pe_Qe(RR0, Stif, Fc, Oe, Pe, Qe) - REAL(BDKi), intent(in) :: RR0(:,:), Stif(:,:), Fc(:) + subroutine Calc_Oe_Pe_Qe(RR0, Stif, E1, Fc, cet, k1s, Oe, Pe, Qe) + REAL(BDKi), intent(in) :: RR0(:,:), Stif(:,:), E1(:), Fc(:), cet, k1s REAL(BDKi), intent(inout) :: Oe(:,:), Pe(:,:), Qe(:,:) REAL(BDKi) :: Wrk33(3,3) REAL(BDKi) :: tildeE(3,3) @@ -2750,7 +2761,7 @@ subroutine Calc_Oe_Pe_Qe(RR0, Stif, Fc, Oe, Pe, Qe) Wrk33 = OuterProduct(RR0(1:3,3), RR0(1:3,3)) ! z-direction in IEC coords C21 = Stif(4:6,1:3) + cet*k1s*Wrk33(:,:) - tildeE = SkewSymMat(m%qp%E1(:,idx_qp,nelem)) + tildeE = SkewSymMat(E1) epsi = MATMUL(Stif(1:3,1:3),tildeE) ! Stif is RR0 * p%Stif0_QP * RR0^T mu = MATMUL(C21,tildeE) @@ -2793,9 +2804,9 @@ subroutine Calc_Oe_Pe_Qe(RR0, Stif, Fc, Oe, Pe, Qe) Qe(4:6,4:6) = -MATMUL(tildeE,Oe(1:3,4:6)) end subroutine - subroutine Calc_Fc_Fd(RR0, uuu, E1, Fc, Fd) - REAL(BDKi), intent(in) :: RR0(:,:), uuu(:), E1(:) - REAL(BDKi), intent(inout) :: Fc(:), Fd(:) + subroutine Calc_Fc_Fd(RR0, uuu, E10, E1, kappa, Stif0, Stif, Fc, Fd, cet, k1s) + REAL(BDKi), intent(in) :: RR0(:,:), uuu(:), E10(:), E1(:), kappa(:), Stif0(:,:), Stif(:,:) + REAL(BDKi), intent(out) :: Fc(:), Fd(:), cet, k1s REAL(BDKi) :: e1s REAL(BDKi) :: eee(6) !< intermediate array for calculation Strain and curvature terms of Fc REAL(BDKi) :: fff(6) !< intermediate array for calculation of the elastic force, Fc @@ -2815,7 +2826,7 @@ subroutine Calc_Fc_Fd(RR0, uuu, E1, Fc, Fd) !! and the transpose for the other direction. ! eee(1:3) = m%qp%E1(1:3,idx_qp,nelem) - m%qp%RR0(1:3,3,idx_qp,nelem) ! Using RR0 z direction in IEC coords call BD_CrvMatrixR(uuu(4:6), R) ! Get rotation at QP as a matrix - Rx0p = matmul(R,p%E10(:,idx_qp,nelem)) ! Calculate rotated initial tangent + Rx0p = matmul(R,E10) ! Calculate rotated initial tangent eee(1:3) = E1(1:3) - Rx0p ! Use rotated initial tangent in place of RR0*i1 to eliminate likely mismatch between R0*i1 and x0' !> ### Set the 1D sectional curvature, \f$ \underline{\kappa} \f$, equation (5) @@ -2836,7 +2847,7 @@ subroutine Calc_Fc_Fd(RR0, uuu, E1, Fc, Fd) !! \f$ !! In other words, \f$ \tilde{k} = \left(\underline{\underline{R}}^\prime\underline{\underline{R}}^T \right) \f$. !! Note: \f$ \underline{\kappa} \f$ was already calculated in the BD_DisplacementQP routine - eee(4:6) = m%qp%kappa(1:3,idx_qp,nelem) + eee(4:6) = kappa(1:3) !FIXME: note that the k_i terms may not be documented correctly here. @@ -2868,7 +2879,7 @@ subroutine Calc_Fc_Fd(RR0, uuu, E1, Fc, Fd) !! \underline{k} !! \end{array} \right\} \f$ !! - fff(1:6) = MATMUL(m%qp%Stif(:,:,idx_qp,nelem),eee) + fff(1:6) = MATMUL(Stif,eee) !> ###Calculate the extension twist coupling. @@ -2903,7 +2914,7 @@ subroutine Calc_Fc_Fd(RR0, uuu, E1, Fc, Fd) !! Note that with coverting to the FAST / IEC coordinate system, we now are using the Ixx and Iyy terms which are located at !! \f$ C_{et} = C_{4,4} + C_{5,5} \f$ ! Refer Section 1.4 in "Dymore User's Manual - Formulation and finite element implementation of beam elements". - cet= p%Stif0_QP(4,4,(nelem-1)*p%nqp+idx_qp) + p%Stif0_QP(5,5,(nelem-1)*p%nqp+idx_qp) ! Dymore theory (22) + cet = Stif0(4,4) + Stif0(5,5) ! Dymore theory (22) Fc(1:3) = fff(1:3) + 0.5_BDKi*cet*k1s*k1s*RR0(1:3,3) ! Dymore theory (25a). Note z-axis is the length of blade. Fc(4:6) = fff(4:6) + cet*e1s*k1s*RR0(1:3,3) ! Dymore theory (25b). Note z-axis is the length of blade. @@ -3085,6 +3096,7 @@ SUBROUTINE BD_DissipativeForce( nelem, p, m,fact ) LOGICAL, INTENT(IN ) :: fact INTEGER(IntKi) :: idx_qp !< index of current quadrature point + REAL(BDKi) :: ffd_t(6) IF (.NOT. fact) then ! skip all but Fc and Fd terms @@ -3095,7 +3107,8 @@ SUBROUTINE BD_DissipativeForce( nelem, p, m,fact ) m%qp%vvp(:,idx_qp,nelem), & m%qp%betaC(:,:,idx_qp,nelem), & m%qp%Fc(:,idx_qp,nelem), & - m%qp%Fd(:,idx_qp,nelem)) + m%qp%Fd(:,idx_qp,nelem), & + ffd_t) END DO ! bjj: we don't use these values when fact is FALSE, so let's save time and ignore them here, too. @@ -3117,11 +3130,13 @@ SUBROUTINE BD_DissipativeForce( nelem, p, m,fact ) m%qp%vvp(:,idx_qp,nelem), & m%qp%betaC(:,:,idx_qp,nelem), & m%qp%Fc(:,idx_qp,nelem), & - m%qp%Fd(:,idx_qp,nelem)) + m%qp%Fd(:,idx_qp,nelem), & + ffd_t) call Calc_Sd_Pd_Od_Qd_Gd_Xd_Yd(m%qp%E1(:,idx_qp,nelem), & m%qp%vvp(:,idx_qp,nelem), & m%qp%betaC(:,:,idx_qp,nelem), & + ffd_t, & m%qp%Sd(:,:,idx_qp,nelem), & m%qp%Od(:,:,idx_qp,nelem), & m%qp%Qd(:,:,idx_qp,nelem), & @@ -3133,14 +3148,13 @@ SUBROUTINE BD_DissipativeForce( nelem, p, m,fact ) ENDIF CONTAINS - subroutine Calc_Sd_Pd_Od_Qd_Gd_Xd_Yd(E1, vvp, betaC, Sd, Od, Qd, Gd, Xd, Yd, Pd) - REAL(BDKi), intent(in) :: E1(:), vvp(:), betaC(:,:) - REAL(BDKi), intent(inout) :: Sd(:,:), Od(:,:), Qd(:,:), Gd(:,:), Xd(:,:), Yd(:,:), Pd(:,:) + subroutine Calc_Sd_Pd_Od_Qd_Gd_Xd_Yd(E1, vvp, betaC, ffd, Sd, Od, Qd, Gd, Xd, Yd, Pd) + REAL(BDKi), intent(in) :: E1(:), vvp(:), betaC(:,:), ffd(:) + REAL(BDKi), intent(out) :: Sd(:,:), Od(:,:), Qd(:,:), Gd(:,:), Xd(:,:), Yd(:,:), Pd(:,:) REAL(BDKi) :: D11(3,3), D12(3,3), D21(3,3), D22(3,3) REAL(BDKi) :: b11(3,3), b12(3,3) REAL(BDKi) :: alpha(3,3) REAL(BDKi) :: SS_ome(3,3) - REAL(BDKi) :: ffd(6) D11 = betaC(1:3,1:3) D12 = betaC(1:3,4:6) @@ -3188,10 +3202,10 @@ subroutine Calc_Sd_Pd_Od_Qd_Gd_Xd_Yd(E1, vvp, betaC, Sd, Od, Qd, Gd, Xd, Yd, Pd) Yd(4:6,4:6) = b12 end subroutine - SUBROUTINE Calc_FC_FD_ffd(E1, vvv, vvp, betaC, Fc, Fd) + SUBROUTINE Calc_FC_FD_ffd(E1, vvv, vvp, betaC, Fc, Fd, ffd) REAL(BDKi), intent(in) :: E1(:), vvv(:), vvp(:), betaC(:,:) - REAL(BDKi), intent(inout) :: Fc(:), Fd(:) - REAL(BDKi) :: eed(6), ffd(6) + REAL(BDKi), intent(out) :: Fc(:), Fd(:), ffd(:) + REAL(BDKi) :: eed(6) ! Compute strain rates eed = vvp From f93ef0504765569fd7c0038ada621d0f64d99b38 Mon Sep 17 00:00:00 2001 From: Derek Slaughter Date: Thu, 29 Aug 2024 18:41:46 +0000 Subject: [PATCH 05/10] Modfy executeOpenfastLinearRegressionCase.py to skip comparison of numbers less than 5e-5 --- reg_tests/executeOpenfastLinearRegressionCase.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/reg_tests/executeOpenfastLinearRegressionCase.py b/reg_tests/executeOpenfastLinearRegressionCase.py index fbdbe635e3..f8c2aa191d 100644 --- a/reg_tests/executeOpenfastLinearRegressionCase.py +++ b/reg_tests/executeOpenfastLinearRegressionCase.py @@ -292,6 +292,9 @@ def indent(msg, sindent='\t'): # Loop through elements where Mloc is not within tolerance of Mbas # Retest to get error message for n, (i,j) in enumerate(zip(*np.where(M_in_tol == False)), 1): + # Skip comparisons of really small numbers + if np.abs(Mloc[i,j]) < 5e-5 and np.abs(Mloc[i,j]) < 5e-5: + continue try: np.testing.assert_allclose(Mloc[i,j], Mbas[i,j], rtol=rtol, atol=atol) except Exception as e: From 63c9d3beff5c341dbf5079622fbcb827283c8c6b Mon Sep 17 00:00:00 2001 From: Derek Slaughter Date: Wed, 18 Sep 2024 13:54:32 +0000 Subject: [PATCH 06/10] Revert "Modfy executeOpenfastLinearRegressionCase.py to skip comparison of numbers less than 5e-5" This reverts commit f93ef0504765569fd7c0038ada621d0f64d99b38. --- reg_tests/executeOpenfastLinearRegressionCase.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/reg_tests/executeOpenfastLinearRegressionCase.py b/reg_tests/executeOpenfastLinearRegressionCase.py index f8c2aa191d..fbdbe635e3 100644 --- a/reg_tests/executeOpenfastLinearRegressionCase.py +++ b/reg_tests/executeOpenfastLinearRegressionCase.py @@ -292,9 +292,6 @@ def indent(msg, sindent='\t'): # Loop through elements where Mloc is not within tolerance of Mbas # Retest to get error message for n, (i,j) in enumerate(zip(*np.where(M_in_tol == False)), 1): - # Skip comparisons of really small numbers - if np.abs(Mloc[i,j]) < 5e-5 and np.abs(Mloc[i,j]) < 5e-5: - continue try: np.testing.assert_allclose(Mloc[i,j], Mbas[i,j], rtol=rtol, atol=atol) except Exception as e: From 56a97f6e4094b228a2f4c69e2e9defd0e2cd7456 Mon Sep 17 00:00:00 2001 From: Derek Slaughter Date: Wed, 18 Sep 2024 14:23:58 +0000 Subject: [PATCH 07/10] Perturb BD rotation states in WM parameters --- modules/beamdyn/src/BeamDyn_IO.f90 | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/modules/beamdyn/src/BeamDyn_IO.f90 b/modules/beamdyn/src/BeamDyn_IO.f90 index 08a1f714fb..15ed92d7d6 100644 --- a/modules/beamdyn/src/BeamDyn_IO.f90 +++ b/modules/beamdyn/src/BeamDyn_IO.f90 @@ -2472,6 +2472,7 @@ SUBROUTINE Perturb_x( p, fieldIndx, node, dof, perturb_sign, x, dx ) REAL(R8Ki) :: orientation(3,3) REAL(R8Ki) :: rotation(3,3) + REAL(R8Ki) :: CrvPerturb(3), CrvBase(3) dx = p%dx(dof) @@ -2479,13 +2480,16 @@ SUBROUTINE Perturb_x( p, fieldIndx, node, dof, perturb_sign, x, dx ) if (dof < 4) then ! translational displacement x%q( dof, node ) = x%q( dof, node ) + dx * perturb_sign else ! w-m parameters - call BD_CrvMatrixR( x%q( 4:6, node ), rotation ) ! returns the rotation matrix (transpose of DCM) that was stored in the state as a w-m parameter - orientation = transpose(rotation) - CALL PerturbOrientationMatrix( orientation, dx * perturb_sign, dof-3 ) ! NOTE: call not using DCM_logmap - - rotation = transpose(orientation) - call BD_CrvExtractCrv( rotation, x%q( 4:6, node ), ErrStat2, ErrMsg2 ) ! return the w-m parameters of the new orientation + ! Calculate perturbation in WM parameters + CrvPerturb = 0.0_R8Ki + CrvPerturb(dof-3) = 4.0_R8Ki * tan(dx * perturb_sign / 4.0_R8Ki) + + ! Get base rotation in WM parameters + CrvBase = x%q(4:6, node) + + ! Compose pertubation and base rotation and store in state + call BD_CrvCompose(x%q(4:6, node), CrvPerturb, CrvBase, FLAG_R1R2) end if else x%dqdt( dof, node ) = x%dqdt( dof, node ) + dx * perturb_sign From 68423c362d5d144a529b7910bbb99cb5be1012e8 Mon Sep 17 00:00:00 2001 From: Derek Slaughter Date: Fri, 20 Sep 2024 15:34:01 +0000 Subject: [PATCH 08/10] BeamDyn: fix formatting in test_tools.F90 --- modules/beamdyn/tests/test_tools.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/modules/beamdyn/tests/test_tools.F90 b/modules/beamdyn/tests/test_tools.F90 index 4850f1770b..3dfa5e31c4 100644 --- a/modules/beamdyn/tests/test_tools.F90 +++ b/modules/beamdyn/tests/test_tools.F90 @@ -173,8 +173,8 @@ type(BD_MiscVarType) function simpleMiscVarType(nqp, dof_node, elem_total, nodes call AllocAry(m%qp%RR0mEta, 3, nqp, elem_total, 'qp_RR0mEta', ErrStat, ErrMsg) call AllocAry(m%DistrLoad_QP, 6, nqp, elem_total, 'DistrLoad_QP', ErrStat, ErrMsg) - call AllocAry(m%qp%uuu, dof_node, nqp, elem_total, 'm%qp%uuu displacement at quadrature point', ErrStat, ErrMsg) - call AllocAry(m%qp%uup, dof_node, nqp, elem_total, 'm%qp%uup displacement prime at quadrature point', ErrStat, ErrMsg) + call AllocAry(m%qp%uuu, dof_node, nqp, elem_total, 'm%qp%uuu displacement at quadrature point', ErrStat, ErrMsg) + call AllocAry(m%qp%uup, dof_node, nqp, elem_total, 'm%qp%uup displacement prime at quadrature point', ErrStat, ErrMsg) ! E1, kappa -- used in force calculations CALL AllocAry(m%qp%E1, dof_node/2,nqp,elem_total, 'm%qp%E1 at quadrature point',ErrStat,ErrMsg) From 11f1cd332bc419ae4296d2684f2653709ed692b1 Mon Sep 17 00:00:00 2001 From: Derek Slaughter Date: Fri, 20 Sep 2024 15:41:59 +0000 Subject: [PATCH 09/10] Update r-test pointer --- reg_tests/r-test | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/reg_tests/r-test b/reg_tests/r-test index 548bf21436..a8ddb774a0 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 548bf2143633a7066f4582167fdf4180951fd7f6 +Subproject commit a8ddb774a0724e0880e95f32b7a9db628aaaf2a8 From e1d3cdb0ce1607504a45c239278cba8bc0fc872c Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Fri, 20 Sep 2024 12:50:19 -0600 Subject: [PATCH 10/10] BD: change to LAPACK_GEMM, add note on ignored error handling --- modules/beamdyn/src/BeamDyn.f90 | 34 ++++++++++++++++++--------------- 1 file changed, 19 insertions(+), 15 deletions(-) diff --git a/modules/beamdyn/src/BeamDyn.f90 b/modules/beamdyn/src/BeamDyn.f90 index 172fca1bd9..25eff3875b 100644 --- a/modules/beamdyn/src/BeamDyn.f90 +++ b/modules/beamdyn/src/BeamDyn.f90 @@ -2360,8 +2360,8 @@ SUBROUTINE BD_DisplacementQP( nelem, p, x, m ) TYPE(BD_ContinuousStateType), INTENT(IN ) :: x !< Continuous states at t TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< misc/optimization variables - INTEGER(IntKi) :: ErrStat !< index to current element - CHARACTER(ErrMsgLen) :: ErrMsg !< index to current element + INTEGER(IntKi) :: ErrStat !< Ignored error handling for LAPACK_GEMM + CHARACTER(ErrMsgLen) :: ErrMsg !< Ignored error handling for LAPACK_GEMM INTEGER(IntKi) :: idx_qp !< index to the current quadrature point INTEGER(IntKi) :: elem_start !< Node point of first node in current element @@ -2387,8 +2387,9 @@ SUBROUTINE BD_DisplacementQP( nelem, p, x, m ) elem_start = p%node_elem_idx(nelem,1) ! Use matrix multiplication to interpolate position and position derivative to quadrature points - call LAPACK_DGEMM('N','N', 1.0_BDKi, x%q(1:3,elem_start:elem_start+p%nodes_per_elem-1), p%Shp, 0.0_BDKi, m%qp%uuu(1:3,:,nelem), ErrStat, ErrMsg) - call LAPACK_DGEMM('N','N', 1.0_BDKi, x%q(1:3,elem_start:elem_start+p%nodes_per_elem-1), p%ShpDer, 0.0_BDKi, m%qp%uup(1:3,:,nelem), ErrStat, ErrMsg) + ! NOTE: errors from LAPACK_GEMM can only be due to matrix size mismatch, so they can be safely ignored if matrices are correct size + call LAPACK_GEMM('N','N', 1.0_BDKi, x%q(1:3,elem_start:elem_start+p%nodes_per_elem-1), p%Shp, 0.0_BDKi, m%qp%uuu(1:3,:,nelem), ErrStat, ErrMsg) + call LAPACK_GEMM('N','N', 1.0_BDKi, x%q(1:3,elem_start:elem_start+p%nodes_per_elem-1), p%ShpDer, 0.0_BDKi, m%qp%uup(1:3,:,nelem), ErrStat, ErrMsg) ! Apply Jacobian to get position derivative with respect to X-axis do idx_qp = 1, p%nqp @@ -2414,8 +2415,8 @@ SUBROUTINE BD_RotationalInterpQP( nelem, p, x, m ) TYPE(BD_ContinuousStateType), INTENT(IN ) :: x !< Continuous states at t TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< misc/optimization variables - INTEGER(IntKi) :: ErrStat !< index to current element - CHARACTER(ErrMsgLen) :: ErrMsg !< index to current element + INTEGER(IntKi) :: ErrStat !< Ignored error handling for LAPACK_GEMM + CHARACTER(ErrMsgLen) :: ErrMsg !< Ignored error handling for LAPACK_GEMM INTEGER(IntKi) :: idx_qp !< index to the current quadrature point INTEGER(IntKi) :: elem_start !< Node point of first node in current element INTEGER(IntKi) :: idx_node !< index to current GLL point in element @@ -2480,8 +2481,9 @@ SUBROUTINE BD_RotationalInterpQP( nelem, p, x, m ) ! Use matrix multiplication to interpolate rotation and rotation derivative to quadrature points ! These rotations do not include the root node rotation at this point (added later in function) - call LAPACK_DGEMM('N','N', 1.0_BDKi, m%Nrrr(:,:,nelem), p%Shp, 0.0_BDKi, m%qp%uuu(4:6,:,nelem), ErrStat, ErrMsg) - call LAPACK_DGEMM('N','N', 1.0_BDKi, m%Nrrr(:,:,nelem), p%ShpDer, 0.0_BDKi, m%qp%uup(4:6,:,nelem), ErrStat, ErrMsg) + ! NOTE: errors from LAPACK_GEMM can only be due to matrix size mismatch, so they can be safely ignored if matrices are correct size + call LAPACK_GEMM('N','N', 1.0_BDKi, m%Nrrr(:,:,nelem), p%Shp, 0.0_BDKi, m%qp%uuu(4:6,:,nelem), ErrStat, ErrMsg) + call LAPACK_GEMM('N','N', 1.0_BDKi, m%Nrrr(:,:,nelem), p%ShpDer, 0.0_BDKi, m%qp%uup(4:6,:,nelem), ErrStat, ErrMsg) ! Apply Jacobian to get rotation derivative with respect to X-axis do idx_qp = 1, p%nqp @@ -2950,8 +2952,8 @@ SUBROUTINE BD_QPDataVelocity( p, x, m ) TYPE(BD_ContinuousStateType), INTENT(IN ) :: x !< Continuous states at t TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< Misc/optimization variables - INTEGER(IntKi) :: ErrStat !< index to current element - CHARACTER(ErrMsgLen) :: ErrMsg !< index to current element + INTEGER(IntKi) :: ErrStat !< Ignored error handling for LAPACK_GEMM + CHARACTER(ErrMsgLen) :: ErrMsg !< Ignored error handling for LAPACK_GEMM INTEGER(IntKi) :: nelem !< index to current element INTEGER(IntKi) :: idx_qp !< index to quadrature point INTEGER(IntKi) :: elem_start !< Starting quadrature point of current element @@ -2965,8 +2967,9 @@ SUBROUTINE BD_QPDataVelocity( p, x, m ) elem_start = p%node_elem_idx(nelem,1) ! Use matrix multiplication to interpolate velocity and velocity derivative to quadrature points - call LAPACK_DGEMM('N','N', 1.0_BDKi, x%dqdt(:,elem_start:elem_start+p%nodes_per_elem-1), p%Shp, 0.0_BDKi, m%qp%vvv(:,:,nelem), ErrStat, ErrMsg) - call LAPACK_DGEMM('N','N', 1.0_BDKi, x%dqdt(:,elem_start:elem_start+p%nodes_per_elem-1), p%ShpDer, 0.0_BDKi, m%qp%vvp(:,:,nelem), ErrStat, ErrMsg) + ! NOTE: errors from LAPACK_GEMM can only be due to matrix size mismatch, so they can be safely ignored if matrices are correct size + call LAPACK_GEMM('N','N', 1.0_BDKi, x%dqdt(:,elem_start:elem_start+p%nodes_per_elem-1), p%Shp, 0.0_BDKi, m%qp%vvv(:,:,nelem), ErrStat, ErrMsg) + call LAPACK_GEMM('N','N', 1.0_BDKi, x%dqdt(:,elem_start:elem_start+p%nodes_per_elem-1), p%ShpDer, 0.0_BDKi, m%qp%vvp(:,:,nelem), ErrStat, ErrMsg) ! Apply Jacobian to get velocity derivative with respect to X-axis do idx_qp = 1, p%nqp @@ -2992,8 +2995,8 @@ SUBROUTINE BD_QPDataAcceleration( p, OtherState, m ) TYPE(BD_OtherStateType), INTENT(IN ) :: OtherState !< Other states at t on input; at t+dt on outputs TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< Misc/optimization variables - INTEGER(IntKi) :: ErrStat !< index to current element - CHARACTER(ErrMsgLen) :: ErrMsg !< index to current element + INTEGER(IntKi) :: ErrStat !< Ignored error handling for LAPACK_GEMM + CHARACTER(ErrMsgLen) :: ErrMsg !< Ignored error handling for LAPACK_GEMM INTEGER(IntKi) :: nelem !< index of current element INTEGER(IntKi) :: idx_qp !< index of current quadrature point INTEGER(IntKi) :: idx_node @@ -3005,7 +3008,8 @@ SUBROUTINE BD_QPDataAcceleration( p, OtherState, m ) elem_start = p%node_elem_idx(nelem,1) ! Interpolate the acceleration term at t+dt (OtherState%acc is at t+dt) to quadrature points - call LAPACK_DGEMM('N','N', 1.0_BDKi, OtherState%acc(:,elem_start:elem_start+p%nodes_per_elem-1), p%Shp, 0.0_BDKi, m%qp%aaa(:,:,nelem), ErrStat, ErrMsg) + ! NOTE: errors from LAPACK_GEMM can only be due to matrix size mismatch, so they can be safely ignored if matrices are correct size + call LAPACK_GEMM('N','N', 1.0_BDKi, OtherState%acc(:,elem_start:elem_start+p%nodes_per_elem-1), p%Shp, 0.0_BDKi, m%qp%aaa(:,:,nelem), ErrStat, ErrMsg) end do