From 75f1e32918a156281f1ff7bcd933d443616e6d73 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Fri, 19 Jul 2024 10:58:09 -0600 Subject: [PATCH 01/22] Fix docs -- fails on rtd --- docs/source/install/index.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/source/install/index.rst b/docs/source/install/index.rst index 9c4a4fa46..ee7951a10 100644 --- a/docs/source/install/index.rst +++ b/docs/source/install/index.rst @@ -208,11 +208,11 @@ To pull a specific release, substitute the version number in place of `latest` i Build your own images --------------------- -You can also build your own custom images using our `Dockerfile` or base your images on ours. See -`here `_ for more information on this. +You can also build your own custom images using our `Dockerfile` or base your images on ours. See the +`Docker readme `_ for more information on this. -.. _python_wrapper +.. _python_wrapper: Install the ``openfast_io`` python wrapper ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -229,7 +229,7 @@ or poetry add openfast_io -For more information and installation options, see `here `_. +For more information and installation options, see the `OpenFAST Python readme `_. .. _compile_from_source: From 0e7774a668dd0dbf9e451402dcf1ad0e37279a33 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Fri, 19 Jul 2024 12:47:56 -0600 Subject: [PATCH 02/22] Add !$ OMP critical around file opening for VTK We were having problems with the `!OMP` directives around high resolution file reading in AWAE.f90. Since the file reading starts with a call to `GetNewUnit` before starting the opening, parallel calls to `GetNewUnit` could result in the same unit number handed out to two processes. The first process would open the file and start reading, but then the second process would open a different file with the same unit number causing read errors for both processes as they attempted to read the same file at the same time. Adding `$OMP critical` around the `GetNewUnit` and following `OpenFile...` calls so that in theory these cannot be done in parallel. Co-authored-by: Derek Slaughter --- modules/nwtc-library/src/VTK.f90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/modules/nwtc-library/src/VTK.f90 b/modules/nwtc-library/src/VTK.f90 index d70345dec..fb57a5781 100644 --- a/modules/nwtc-library/src/VTK.f90 +++ b/modules/nwtc-library/src/VTK.f90 @@ -158,8 +158,10 @@ SUBROUTINE ReadVTK_SP_info( FileName, descr, dims, origin, gridSpacing, vecLabel closeOnReturn = .FALSE. END IF + !$OMP critical CALL GetNewUnit( Un, ErrStat, ErrMsg ) CALL OpenFInpFile ( Un, TRIM(FileName), ErrStat, ErrMsg ) + !$OMP end critical if (ErrStat >= AbortErrLev) return CALL ReadCom( Un, FileName, 'File header: Module Version (line 1)', ErrStat2, ErrMsg2, 0 ) @@ -358,8 +360,10 @@ SUBROUTINE WrVTK_SP_header( FileName, descr, Un, ErrStat, ErrMsg ) INTEGER(IntKi) , INTENT( OUT) :: ErrStat !< error level/status of OpenFOutFile operation CHARACTER(*) , INTENT( OUT) :: ErrMsg !< message when error occurs + !$OMP critical CALL GetNewUnit( Un, ErrStat, ErrMsg ) CALL OpenFOutFile ( Un, TRIM(FileName), ErrStat, ErrMsg ) + !$OMP end critical if (ErrStat >= AbortErrLev) return WRITE(Un,'(A)') '# vtk DataFile Version 3.0' From 97e0f5bddcb602d51be1fbc2e9d4ff528436ba62 Mon Sep 17 00:00:00 2001 From: Guillaume Jacquenot Date: Sat, 10 Aug 2024 22:22:04 +0200 Subject: [PATCH 03/22] :pencil: Fixed typo in docs/source/install/index.rst --- docs/source/install/index.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/install/index.rst b/docs/source/install/index.rst index 85c33ac20..dfab67c3a 100644 --- a/docs/source/install/index.rst +++ b/docs/source/install/index.rst @@ -159,7 +159,7 @@ containing the executables, and running a simple test command: Running OpenFAST with docker ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -OpenFAST is avilable to be run on docker starting with version 3.5.3. Three approaches are shared below. +OpenFAST is available to be run on docker starting with version 3.5.3. Three approaches are shared below. Using a docker image from Docker hub ------------------------------------ From ff67ad42966eb20ad95a753c6192c20950a2f27d Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Fri, 23 Aug 2024 14:11:27 -0600 Subject: [PATCH 04/22] FileInfoType: increase line length allowed The old line length was 1024 characters. Increased this to 8192 -- May be need in HD if a large number of floating bodies are specified --- modules/nwtc-library/src/NWTC_Base.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modules/nwtc-library/src/NWTC_Base.f90 b/modules/nwtc-library/src/NWTC_Base.f90 index 2ccec45fb..b27652c69 100644 --- a/modules/nwtc-library/src/NWTC_Base.f90 +++ b/modules/nwtc-library/src/NWTC_Base.f90 @@ -40,7 +40,7 @@ MODULE NWTC_Base INTEGER(IntKi), PARAMETER :: MinChanLen = 10 !< The min allowable length of channel names (i.e., width of output columns), used because some modules (like Bladed DLL outputs) have excessively long names INTEGER(IntKi), PARAMETER :: LinChanLen = 200 !< The allowable length of row/column names in linearization files - INTEGER(IntKi), PARAMETER :: MaxFileInfoLineLen = 1024 !< The allowable length of an input line stored in FileInfoType%Lines + INTEGER(IntKi), PARAMETER :: MaxFileInfoLineLen = 8192 !< The allowable length of an input line stored in FileInfoType%Lines INTEGER(IntKi), PARAMETER :: NWTC_Verbose = 10 !< The maximum level of verbosity INTEGER(IntKi), PARAMETER :: NWTC_VerboseLevel = 5 !< a number in [0, NWTC_Verbose]: 0 = no output; NWTC_Verbose=verbose; From d97ee668c5ca3f83827485ab50f9370ad3e590fd Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Wed, 28 Aug 2024 12:57:16 -0600 Subject: [PATCH 05/22] HD bug fix: prevent array index out-of-bound error in HDOut_MapOutputs when more than 9 potential-flow bodies are present --- modules/hydrodyn/src/HydroDyn_Output.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/modules/hydrodyn/src/HydroDyn_Output.f90 b/modules/hydrodyn/src/HydroDyn_Output.f90 index 44668a17c..a48f6429e 100644 --- a/modules/hydrodyn/src/HydroDyn_Output.f90 +++ b/modules/hydrodyn/src/HydroDyn_Output.f90 @@ -1161,7 +1161,7 @@ SUBROUTINE HDOut_MapOutputs( CurrentTime, p, y, m_WAMIT, m_WAMIT2, NWaveElev, Wa ! Need to use individual components of force for output reporting, the y%mesh data has total forces from all contributions if ( p%potMod == 1 ) then if ( p%NBodyMod == 1 .or. p%NBody == 1 ) then - do iBody = 1,p%NBody + do iBody = 1,min(p%NBody,9) ! Can only output the first 9 bodies for now startIndx = 6*(iBody-1) + 1 endIndx = startIndx + 5 AllOuts(FWaves1 (:,iBody)) = m_WAMIT(1)%F_Waves1(startIndx:endIndx) @@ -1179,7 +1179,7 @@ SUBROUTINE HDOut_MapOutputs( CurrentTime, p, y, m_WAMIT, m_WAMIT2, NWaveElev, Wa ! This happens when NBodyMod > 1, in which case, each WAMIT object is for a single body, but there may be multiple bodies in the HydroDyn model, ! so we need to use BodyID to determine the index into the complete HydroDyn list of WAMIT bodies - do iBody = 1,p%NBody + do iBody = 1,min(p%NBody,9) ! Can only output the first 9 bodies for now startIndx = 6*(iBody-1) + 1 endIndx = startIndx + 5 AllOuts(FWaves1 (:,iBody)) = m_WAMIT(iBody)%F_Waves1 From 6253e3e99ec91a179de20e344117dee789bb901f Mon Sep 17 00:00:00 2001 From: Derek Slaughter Date: Fri, 16 Aug 2024 14:36:32 +0000 Subject: [PATCH 06/22] 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 2f8dbfd1e..ad69f8747 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 519a40e58..c1b9796a0 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 666137035..7ede50f62 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 df2723a60..6d9a7325d 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 2daf088cc..8d7eb9338 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 07/22] 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 ad69f8747..d37909b0e 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 08/22] 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 1f64ec584..4850f1770 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 09/22] 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 d37909b0e..172fca1bd 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 10/22] 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 fbdbe635e..f8c2aa191 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 11/22] 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 f8c2aa191..fbdbe635e 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 12/22] 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 08a1f714f..15ed92d7d 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 13/22] 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 4850f1770..3dfa5e31c 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 14/22] 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 548bf2143..a8ddb774a 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 15/22] 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 172fca1bd..25eff3875 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 From dde75b8a388829395c5b4b44dc23dcf7ed78fc4e Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Mon, 7 Oct 2024 18:32:19 -0600 Subject: [PATCH 16/22] Intel builds: add /heap-arrays:1000 compile flag We have been having issues with stack overflows in FAST.Farm when large wind grids were passed to AWAE. This was eventually tracked down to line 1078 in AWAE.f90 which reads: `m%u_IfW_Low%PositionXYZ = p%Grid_low` The `p%Grid_low` is of unknown size at compile time, but can be extremely large (3x160000 or more). This copy involves a temporary array and would normally be handled on the stack, but could result in an overflow for some models. By setting the `/heap-arrays:1000` any operation resulting in a temporary array of unknown array size at compile time will use the heap for the temprary array, and as will any array known at compile time to be larger than 1000 kB. Testing shows that this fixes issue #2053, and will likely also solve #843 and #2241 See https://www.intel.com/content/www/us/en/docs/fortran-compiler/developer-guide-reference/2024-2/heap-arrays.html for reference. --- vs-build/FAST-farm/FAST-Farm.vfproj | 6 +++--- vs-build/FASTlib/FASTlib.vfproj | 14 +++++++------- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/vs-build/FAST-farm/FAST-Farm.vfproj b/vs-build/FAST-farm/FAST-Farm.vfproj index e7c3152f3..5f1d3153b 100644 --- a/vs-build/FAST-farm/FAST-Farm.vfproj +++ b/vs-build/FAST-farm/FAST-Farm.vfproj @@ -25,7 +25,7 @@ - + @@ -35,7 +35,7 @@ - + @@ -55,7 +55,7 @@ - + diff --git a/vs-build/FASTlib/FASTlib.vfproj b/vs-build/FASTlib/FASTlib.vfproj index 7a94c7ae6..32c624b94 100644 --- a/vs-build/FASTlib/FASTlib.vfproj +++ b/vs-build/FASTlib/FASTlib.vfproj @@ -23,7 +23,7 @@ - + @@ -32,7 +32,7 @@ - + @@ -50,7 +50,7 @@ - + @@ -68,7 +68,7 @@ - + @@ -86,7 +86,7 @@ - + @@ -104,7 +104,7 @@ - + @@ -122,7 +122,7 @@ - + From 31d281e444f704d87e1b959efec8b45360c83d57 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Mon, 7 Oct 2024 18:45:41 -0600 Subject: [PATCH 17/22] Add /heap-arrays:1000 to cmake intel builds --- cmake/OpenfastFortranOptions.cmake | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/cmake/OpenfastFortranOptions.cmake b/cmake/OpenfastFortranOptions.cmake index 0efa0aed5..eba55e8df 100644 --- a/cmake/OpenfastFortranOptions.cmake +++ b/cmake/OpenfastFortranOptions.cmake @@ -158,7 +158,9 @@ endmacro(set_fast_intel_fortran) # arch # macro(set_fast_intel_fortran_posix) - set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fpic -fpp") + # Set size where temporary are stored on heap instead of stack + # 1000: size in kB (1 MB) + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fpic -fpp -heap-arrays 1000") # debug flags if(CMAKE_BUILD_TYPE MATCHES Debug) @@ -201,7 +203,9 @@ macro(set_fast_intel_fortran_windows) # Turn off specific warnings # - 5199: too many continuation lines # - 5268: 132 column limit - set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} /Qdiag-disable:5199,5268 /fpp") + # Set size where temporary are stored on heap instead of stack + # 1000: size in kB (1 MB) + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} /Qdiag-disable:5199,5268 /fpp /heap-arrays:1000") # If double precision, make constants double precision if (DOUBLE_PRECISION) From a2e0381bf1abdaf78a856109602a25a5711abb0d Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Mon, 21 Oct 2024 16:28:29 -0600 Subject: [PATCH 18/22] Add changelog for 3.5.4 --- docs/changelogs/v3.5.4.md | 91 +++++++++++++++++++++++++++++++++ docs/source/user/api_change.rst | 6 +++ 2 files changed, 97 insertions(+) create mode 100644 docs/changelogs/v3.5.4.md diff --git a/docs/changelogs/v3.5.4.md b/docs/changelogs/v3.5.4.md new file mode 100644 index 000000000..2cd5054ca --- /dev/null +++ b/docs/changelogs/v3.5.4.md @@ -0,0 +1,91 @@ +**Feature or improvement description** +Pull request to merge `rc-3.5.4` into `main` and create a tagged release for v3.5.4. + +See the milestone and project pages for additional information + + https://github.com/OpenFAST/openfast/milestone/14 + +Test results, if applicable +See GitHub Actions + +### Release checklist: +- [ ] Update the documentation version in docs/conf.py +- [ ] Update the versions in docs/source/user/api_change.rst +- [ ] Verify readthedocs builds correctly +- [ ] Create a tag in OpenFAST +- [ ] Create a merge commit in r-test and add a corresponding annotated tag +- [ ] Compile executables for Windows builds + - [ ] AeroDyn_Driver_x64.exe + - [ ] AeroDyn_Driver_x64_OpenMP.exe + - [ ] AeroDyn_Inflow_C_Binding_x64.dll + - [ ] AeroDyn_Inflow_C_Binding_x64_OpenMP.dll + - [ ] BeamDyn_Driver_x64.exe + - [ ] DISCON.dll (x64) + - [ ] DISCON_ITIBarge.dll (x64) + - [ ] DISCON_OC3Hywind.dll (x64) + - [ ] DISCON_SC.dll (x64) + - [ ] FAST.Farm_x64.exe + - [ ] FAST.Farm_x64_OMP.exe + - [ ] FAST_SFunc.mexw64 + - [ ] HydroDynDriver_x64.exe + - [ ] HydroDyn_C_Binding_x64.dll + - [ ] IfW_C_Binding_x64.dll + - [ ] InflowWind_Driver_x64.exe + - [ ] InflowWind_Driver_x64_OpenMP.exe + - [ ] MoorDyn_Driver_x64.exe + - [ ] MoorDyn_C_Binding_x64.dll + - [ ] OpenFAST-Simulink_x64.dll + - [ ] openfast_x64.exe + - [ ] Turbsim_x64.exe + +# Changelog + +## Overview + +This release includes performance improvements for BeamDyn (up to 30% speed increase), and a fix for stack overflows with FAST.Farm (Intel compilation). A few other minor pdates are inluded as outlined below. + +Anyone using BeamDyn will want to update to this version. + + +## General + +### Build systems + +#2311 Always build `openfastcpplib` as shared. Use `BUILD_OPENFAST_CPP_DRIVER` to disable `openfastcpp` executable (@deslaughter) + +### Docker builds + + + + +## Solvers + + +### FAST.Farm + +#2452 Fix for some stack overflow issues with FAST.Farm when copying large amounts of wind data (closes #2053) (@andrew-platt) +#2340 Add `!$ OMP critical` around file opening for VTK to prevent file collision (@andrew-platt) + +## Module changes + + +### BeamDyn + +#2399 BeamDyn performance improvements (@deslaughter) + + +### ElastoDyn + +#2321 Backport of #2317: Explicitly initialize `ED` `RtHS` to zero -- corrects issue with inccorect linearization results (@andrew-platt) + + + + +## Input file changes + +No input files change with this release, as this only includes minor bugfixes. + +Full list of changes: https://openfast.readthedocs.io/en/main/source/user/api_change.html + +Full input file sets: https://github.com/OpenFAST/r-test/tree/v3.5.4 (example input files from the regression testing) + diff --git a/docs/source/user/api_change.rst b/docs/source/user/api_change.rst index 07af23f90..cf16f2fe4 100644 --- a/docs/source/user/api_change.rst +++ b/docs/source/user/api_change.rst @@ -9,6 +9,12 @@ The changes are tabulated according to the module input file, line number, and f The line number corresponds to the resulting line number after all changes are implemented. Thus, be sure to implement each in order so that subsequent line numbers are correct. +OpenFAST v3.5.3 to OpenFAST v3.5.4 +---------------------------------- + +No input file changes were made. + + OpenFAST v3.5.2 to OpenFAST v3.5.3 ---------------------------------- From 05a16b2672cc1e0b1eb7c7d1f8f68575487583cf Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Mon, 21 Oct 2024 16:43:23 -0600 Subject: [PATCH 19/22] Update version info in conf.py --- docs/conf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/conf.py b/docs/conf.py index 146fb2b41..f1b3d9f22 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -130,7 +130,7 @@ def runDoxygen(sourcfile, doxyfileIn, doxyfileOut): # The short X.Y version. version = u'3.5' # The full version, including alpha/beta/rc tags. -release = u'v3.5.3' +release = u'v3.5.4' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. From ea167509c08e2f559a1458bdf901a618540aaf45 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Mon, 21 Oct 2024 17:17:37 -0600 Subject: [PATCH 20/22] Correction to changelog v3.5.4.md --- docs/changelogs/v3.5.4.md | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/docs/changelogs/v3.5.4.md b/docs/changelogs/v3.5.4.md index 2cd5054ca..8862d40a8 100644 --- a/docs/changelogs/v3.5.4.md +++ b/docs/changelogs/v3.5.4.md @@ -42,7 +42,7 @@ See GitHub Actions ## Overview -This release includes performance improvements for BeamDyn (up to 30% speed increase), and a fix for stack overflows with FAST.Farm (Intel compilation). A few other minor pdates are inluded as outlined below. +This release includes performance improvements for BeamDyn (up to 30% speed increase), python file readers and writers from WEIS, and a fix for stack overflows with FAST.Farm (Intel compilation). A few other minor updates are inluded as outlined below. Anyone using BeamDyn will want to update to this version. @@ -52,22 +52,23 @@ Anyone using BeamDyn will want to update to this version. ### Build systems #2311 Always build `openfastcpplib` as shared. Use `BUILD_OPENFAST_CPP_DRIVER` to disable `openfastcpp` executable (@deslaughter) +#2173 Fix crash in `MAP_End` when using Intel's new icx compiler and disable caching in setup-python GH action (@deslaughter) -### Docker builds +### Python file readers/writers +#2188 Add WEIS file readers and writers (@cortadocodes and @mayankchetan) ## Solvers - ### FAST.Farm #2452 Fix for some stack overflow issues with FAST.Farm when copying large amounts of wind data (closes #2053) (@andrew-platt) #2340 Add `!$ OMP critical` around file opening for VTK to prevent file collision (@andrew-platt) -## Module changes +## Module changes ### BeamDyn From efad7a31e76579f77f44768eb9d5195032a408af Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Thu, 24 Oct 2024 15:13:58 -0600 Subject: [PATCH 21/22] Missing PRs and Typo in v3.5.4.md Thanks to @pablo-benito for noticing the typo --- docs/changelogs/v3.5.4.md | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/docs/changelogs/v3.5.4.md b/docs/changelogs/v3.5.4.md index 8862d40a8..e099bfab5 100644 --- a/docs/changelogs/v3.5.4.md +++ b/docs/changelogs/v3.5.4.md @@ -42,7 +42,7 @@ See GitHub Actions ## Overview -This release includes performance improvements for BeamDyn (up to 30% speed increase), python file readers and writers from WEIS, and a fix for stack overflows with FAST.Farm (Intel compilation). A few other minor updates are inluded as outlined below. +This release includes performance improvements for BeamDyn (up to 30% speed increase), python file readers and writers from WEIS, and a fix for stack overflows with FAST.Farm (Intel compilation). A few other minor updates are included as outlined below. Anyone using BeamDyn will want to update to this version. @@ -80,7 +80,11 @@ Anyone using BeamDyn will want to update to this version. #2321 Backport of #2317: Explicitly initialize `ED` `RtHS` to zero -- corrects issue with inccorect linearization results (@andrew-platt) +### HydroDyn +#2397 HD bug fix: prevent array index out-of-bound error in `HDOut_MapOutputs` when more than 9 potential-flow bodies are present (@luwang00) +### NWTC-Library +#2389 FileInfoType: increase line length allowed (@andrew-platt) ## Input file changes From 8136967120f9fc09ded62b5feb3cb89f64cf4fe7 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Mon, 28 Oct 2024 15:25:22 -0600 Subject: [PATCH 22/22] Update r-test pointer for 3.5.4 --- 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 a8ddb774a..f282c79af 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit a8ddb774a0724e0880e95f32b7a9db628aaaf2a8 +Subproject commit f282c79af7152d20ec7becb7a7c9edbf65cca71c