From 73173d0ba89b05e6f14d54c19a930066eb67f4d0 Mon Sep 17 00:00:00 2001 From: Derek Slaughter Date: Tue, 30 Dec 2025 15:27:59 -0500 Subject: [PATCH 01/12] Remove unused error handling from DCM_logmap --- modules/nwtc-library/src/ModMesh.f90 | 30 ++++---------------- modules/nwtc-library/src/ModMesh_Mapping.f90 | 12 +++----- modules/nwtc-library/src/NWTC_Num.f90 | 22 ++------------ modules/wakedynamics/src/WakeDynamics.f90 | 4 +-- 4 files changed, 14 insertions(+), 54 deletions(-) diff --git a/modules/nwtc-library/src/ModMesh.f90 b/modules/nwtc-library/src/ModMesh.f90 index 4844a765fb..5a9c834c0d 100644 --- a/modules/nwtc-library/src/ModMesh.f90 +++ b/modules/nwtc-library/src/ModMesh.f90 @@ -3205,18 +3205,10 @@ SUBROUTINE MeshExtrapInterp1(u1, u2, tin, u_out, tin_out, ErrStat, ErrMsg ) DO node=1,u_out%Nnodes Orient = u1%Orientation(:,:,node) - CALL DCM_logmap ( Orient, tensor(:,1), ErrStat, ErrMsg ) - IF (ErrStat >= AbortErrLev ) THEN - ErrMsg = 'MeshExtrapInterp1:'//TRIM(ErrMsg) - RETURN - END IF + CALL DCM_logmap(Orient, tensor(:,1)) Orient = u2%Orientation(:,:,node) - CALL DCM_logmap ( Orient, tensor(:,2), ErrStat, ErrMsg ) - IF (ErrStat >= AbortErrLev ) THEN - ErrMsg = 'MeshExtrapInterp1:'//TRIM(ErrMsg) - RETURN - END IF + CALL DCM_logmap(Orient, tensor(:,2)) CALL DCM_SetLogMapForInterp( tensor ) @@ -3347,25 +3339,13 @@ SUBROUTINE MeshExtrapInterp2(u1, u2, u3, tin, u_out, tin_out, ErrStat, ErrMsg ) DO node=1,u_out%Nnodes Orient = u1%Orientation(:,:,node) - CALL DCM_logmap ( Orient, tensor(:,1), ErrStat, ErrMsg ) - IF (ErrStat >= AbortErrLev ) THEN - ErrMsg = 'MeshExtrapInterp2:'//TRIM(ErrMsg) - RETURN - END IF + CALL DCM_logmap(Orient, tensor(:,1)) Orient = u2%Orientation(:,:,node) - CALL DCM_logmap ( Orient, tensor(:,2), ErrStat, ErrMsg ) - IF (ErrStat >= AbortErrLev ) THEN - ErrMsg = 'MeshExtrapInterp2:'//TRIM(ErrMsg) - RETURN - END IF + CALL DCM_logmap(Orient, tensor(:,2)) Orient = u3%Orientation(:,:,node) - CALL DCM_logmap ( Orient, tensor(:,3), ErrStat, ErrMsg ) - IF (ErrStat >= AbortErrLev ) THEN - ErrMsg = 'MeshExtrapInterp2:'//TRIM(ErrMsg) - RETURN - END IF + CALL DCM_logmap(Orient, tensor(:,3)) CALL DCM_SetLogMapForInterp( tensor ) diff --git a/modules/nwtc-library/src/ModMesh_Mapping.f90 b/modules/nwtc-library/src/ModMesh_Mapping.f90 index fd63f374f3..79f3fed31d 100644 --- a/modules/nwtc-library/src/ModMesh_Mapping.f90 +++ b/modules/nwtc-library/src/ModMesh_Mapping.f90 @@ -1333,15 +1333,13 @@ SUBROUTINE Transfer_Motions_Line2_to_Point( Src, Dest, MeshMap, ErrStat, ErrMsg RotationMatrixD = MATMUL( TRANSPOSE( Src%RefOrientation(:,:,n1) ), Src%Orientation(:,:,n1) ) RotationMatrixD = MATMUL( Dest%RefOrientation(:,:,i), RotationMatrixD ) - CALL DCM_logmap( RotationMatrixD, FieldValue(:,1), ErrStat, ErrMsg ) - IF (ErrStat >= AbortErrLev) RETURN + CALL DCM_logmap(RotationMatrixD, FieldValue(:,1)) ! calculate Rotation matrix for FieldValueN2 and convert to tensor: RotationMatrixD = MATMUL( TRANSPOSE( Src%RefOrientation(:,:,n2) ), Src%Orientation(:,:,n2) ) RotationMatrixD = MATMUL( Dest%RefOrientation(:,:,i), RotationMatrixD ) - CALL DCM_logmap( RotationMatrixD, FieldValue(:,2), ErrStat, ErrMsg ) - IF (ErrStat >= AbortErrLev) RETURN + CALL DCM_logmap(RotationMatrixD, FieldValue(:,2)) CALL DCM_SetLogMapForInterp( FieldValue ) ! make sure we don't cross a 2pi boundary @@ -4415,12 +4413,10 @@ SUBROUTINE Create_Augmented_Ln2_Src_Mesh(Src, Dest, MeshMap, Dest_TYPE, ErrStat, ! convert DCMs to tensors: RefOrientationD = Src%RefOrientation(:, :, n1) - CALL DCM_logmap( RefOrientationD, FieldValue(:,1), ErrStat, ErrMsg ) - IF (ErrStat >= AbortErrLev) RETURN + CALL DCM_logmap(RefOrientationD, FieldValue(:,1)) RefOrientationD = Src%RefOrientation(:, :, n2) - CALL DCM_logmap( RefOrientationD, FieldValue(:,2), ErrStat, ErrMsg ) - IF (ErrStat >= AbortErrLev) RETURN + CALL DCM_logmap(RefOrientationD, FieldValue(:,2)) CALL DCM_SetLogMapForInterp( FieldValue ) ! make sure we don't cross a 2pi boundary diff --git a/modules/nwtc-library/src/NWTC_Num.f90 b/modules/nwtc-library/src/NWTC_Num.f90 index d4f36542df..3b73b1865a 100644 --- a/modules/nwtc-library/src/NWTC_Num.f90 +++ b/modules/nwtc-library/src/NWTC_Num.f90 @@ -1258,26 +1258,18 @@ END FUNCTION DCM_expR !! !! This routine is the inverse of DCM_exp (nwtc_num::dcm_exp). \n !! Use DCM_logMap (nwtc_num::dcm_logmap) instead of directly calling a specific routine in the generic interface. - SUBROUTINE DCM_logMapD(DCM, logMap, ErrStat, ErrMsg, thetaOut) + SUBROUTINE DCM_logMapD(DCM, logMap, thetaOut) REAL(R8Ki), INTENT(IN) :: DCM(3,3) !< the direction cosine matrix, \f$\Lambda\f$ REAL(R8Ki), INTENT( OUT) :: logMap(3) !< vector containing \f$\lambda_1\f$, \f$\lambda_2\f$, and \f$\lambda_3\f$, the unique components of skew-symmetric matrix \f$\lambda\f$ REAL(R8Ki),OPTIONAL,INTENT( OUT) :: thetaOut !< the angle of rotation, \f$\theta\f$; output only for debugging - INTEGER(IntKi), INTENT( OUT) :: ErrStat !< Error status of the operation - CHARACTER(*), INTENT( OUT) :: ErrMsg !< Error message if ErrStat /= ErrID_None - ! local variables REAL(R8Ki) :: theta REAL(R8Ki) :: cosTheta REAL(R8Ki) :: TwoSinTheta REAL(R8Ki) :: v(3) REAL(R8Ki) :: divisor - INTEGER(IntKi) :: indx_max - - ! initialization - ErrStat = ErrID_None - ErrMsg = "" - + INTEGER(IntKi) :: indx_max cosTheta = 0.5_DbKi*( trace(DCM) - 1.0_R8Ki ) cosTheta = min( max(cosTheta,-1.0_R8Ki), 1.0_R8Ki ) !make sure it's in a valid range (to avoid cases where this is slightly outside the +/-1 range) @@ -1372,28 +1364,20 @@ SUBROUTINE DCM_logMapD(DCM, logMap, ErrStat, ErrMsg, thetaOut) END SUBROUTINE DCM_logMapD !======================================================================= !> \copydoc nwtc_num::dcm_logmapd - SUBROUTINE DCM_logMapR(DCM, logMap, ErrStat, ErrMsg, thetaOut) + SUBROUTINE DCM_logMapR(DCM, logMap, thetaOut) ! This function computes the logarithmic map for a direction cosine matrix. REAL(SiKi), INTENT(IN) :: DCM(3,3) REAL(SiKi), INTENT( OUT) :: logMap(3) REAL(SiKi),OPTIONAL,INTENT( OUT) :: thetaOut - INTEGER(IntKi), INTENT( OUT) :: ErrStat ! Error status of the operation - CHARACTER(*), INTENT( OUT) :: ErrMsg ! Error message if ErrStat /= ErrID_None - ! local variables REAL(SiKi) :: cosTheta REAL(SiKi) :: theta REAL(SiKi) :: TwoSinTheta REAL(SiKi) :: v(3) REAL(SiKi) :: divisor INTEGER(IntKi) :: indx_max - - ! initialization - ErrStat = ErrID_None - ErrMsg = "" - cosTheta = 0.5_SiKi*( trace(DCM) - 1.0_SiKi ) cosTheta = min( max(cosTheta,-1.0_SiKi), 1.0_SiKi ) !make sure it's in a valid range (to avoid cases where this is slightly outside the +/-1 range) diff --git a/modules/wakedynamics/src/WakeDynamics.f90 b/modules/wakedynamics/src/WakeDynamics.f90 index 8ec6385b4d..dcf119c18a 100644 --- a/modules/wakedynamics/src/WakeDynamics.f90 +++ b/modules/wakedynamics/src/WakeDynamics.f90 @@ -1326,8 +1326,8 @@ subroutine filter_angles2(psi_filt, chi_filt, psi, chi, alpha, alpha_bar) DCM1 = EulerConstruct( (/ psi_filt, 0.0_ReKi, chi_filt /) ) DCM2 = EulerConstruct( (/ psi, 0.0_ReKi, chi /) ) ! Compute the logarithmic map of the DCMs: - CALL DCM_logMap( DCM1, lambda(:,1), errStat, errMsg) - CALL DCM_logMap( DCM2, lambda(:,2), errStat, errMsg) + CALL DCM_logMap(DCM1, lambda(:,1)) + CALL DCM_logMap(DCM2, lambda(:,2)) !Make sure we don't cross a 2pi boundary: CALL DCM_SetLogMapForInterp( lambda ) !Interpolate the logarithmic map: From 92441db9ef318386b68a69f41bc588f326ed5561 Mon Sep 17 00:00:00 2001 From: Derek Slaughter Date: Tue, 30 Dec 2025 15:29:40 -0500 Subject: [PATCH 02/12] Add optional flag to Morison_CalcOutput to skip expensive hydrostatic load calculation. Hydrodyn doesn't need hydrostatic loads when calculating dYdu and perturbing velocities/accelerations. Resulted in 20% speedup. --- modules/hydrodyn/src/HydroDyn.f90 | 27 ++- modules/hydrodyn/src/Morison.f90 | 286 +++++++++++++++++------------- 2 files changed, 185 insertions(+), 128 deletions(-) diff --git a/modules/hydrodyn/src/HydroDyn.f90 b/modules/hydrodyn/src/HydroDyn.f90 index e7825d4264..ff75e252af 100644 --- a/modules/hydrodyn/src/HydroDyn.f90 +++ b/modules/hydrodyn/src/HydroDyn.f90 @@ -1331,7 +1331,7 @@ END SUBROUTINE HydroDyn_UpdateStates !---------------------------------------------------------------------------------------------------------------------------------- !> Routine for computing outputs, used in both loose and tight coupling. -SUBROUTINE HydroDyn_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg ) +SUBROUTINE HydroDyn_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, calcMorisonHstLds ) REAL(DbKi), INTENT(IN ) :: Time !< Current simulation time in seconds TYPE(HydroDyn_InputType), INTENT(INOUT) :: u !< Inputs at Time (note that this is intent out because we're copying the u%WAMITMesh into m%u_wamit%mesh) @@ -1345,6 +1345,8 @@ SUBROUTINE HydroDyn_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, ErrStat, TYPE(HydroDyn_MiscVarType), INTENT(INOUT) :: m !< Initial misc/optimization variables INTEGER(IntKi), INTENT( OUT) :: ErrStat !! Error status of the operation CHARACTER(*), INTENT( OUT) :: ErrMsg !! Error message if ErrStat /= ErrID_None + LOGICAL, OPTIONAL, INTENT(IN ) :: calcMorisonHstLds !< Flag to calculate the Morison hydrostatic loads (default: .true.) + !! Used to speed up Jacobian calculations when perturbing velocity/acceleration inputs INTEGER :: I, J ! Generic counters @@ -1367,11 +1369,18 @@ SUBROUTINE HydroDyn_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, ErrStat, CHARACTER(*), PARAMETER :: RoutineName = 'HydroDyn_CalcOutput' REAL(ReKi), PARAMETER :: LrgAngle = 0.261799387799149 ! Threshold for platform roll and pitch rotation (15 deg). This is consistent with the ElastoDyn check. LOGICAL, SAVE :: FrstWarn_LrgY = .TRUE. + logical :: calcMorisonHstLdsLocal ! Initialize ErrStat ErrStat = ErrID_None ErrMsg = "" + + if (present(calcMorisonHstLds)) then + calcMorisonHstLdsLocal = calcMorisonHstLds + else + calcMorisonHstLdsLocal = .true. + end if ! Write the Hydrodyn-level output file data FROM THE LAST COMPLETED TIME STEP if the user requested module-level output @@ -1567,7 +1576,8 @@ SUBROUTINE HydroDyn_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, ErrStat, IF ( u%Morison%Mesh%Committed ) THEN ! Make sure we are using Morison / there is a valid mesh u%Morison%PtfmRefY = PtfmRefY CALL Morison_CalcOutput( Time, u%Morison, p%Morison, x%Morison, xd%Morison, & - z%Morison, OtherState%Morison, y%Morison, m%Morison, ErrStat2, ErrMsg2 ) + z%Morison, OtherState%Morison, y%Morison, m%Morison, & + ErrStat2, ErrMsg2, calcMorisonHstLdsLocal ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) END IF @@ -1774,6 +1784,7 @@ SUBROUTINE HD_JacobianPInput(Vars, t, u, p, x, xd, z, OtherState, y, m, ErrStat, INTEGER(IntKi) :: i, j, k, col INTEGER(IntKi) :: startingI, startingJ, bOffset, offsetI integer(IntKi) :: iVarWaveElev0, iVarHWindSpeed, iVarPLexp, iVarPropagationDir + logical :: calcMorisonHstLds ErrStat = ErrID_None ErrMsg = '' @@ -1816,19 +1827,27 @@ SUBROUTINE HD_JacobianPInput(Vars, t, u, p, x, xd, z, OtherState, y, m, ErrStat, ! If variable is extended input, skip if (MV_HasFlagsAll(Vars%u(i), VF_ExtLin)) cycle + ! Calculate Morison hydrostatic loads when perturbing displacement/orientation inputs + select case (Vars%u(i)%Field) + case (FieldTransDisp, FieldOrientation) + calcMorisonHstLds = .true. + case default + calcMorisonHstLds = .false. + end select + ! Loop through number of linearization perturbations in variable do j = 1, Vars%u(i)%Num ! Calculate positive perturbation call MV_Perturb(Vars%u(i), j, 1, m%Jac%u, m%Jac%u_perturb) call HydroDyn_VarsUnpackInput(Vars, m%Jac%u_perturb, m%u_perturb) - call HydroDyn_CalcOutput(t, m%u_perturb, p, x, xd, z, OtherState, m%y_lin, m, ErrStat2, ErrMsg2); if (Failed()) return + call HydroDyn_CalcOutput(t, m%u_perturb, p, x, xd, z, OtherState, m%y_lin, m, ErrStat2, ErrMsg2, calcMorisonHstLds); if (Failed()) return call HydroDyn_VarsPackOutput(Vars, m%y_lin, m%Jac%y_pos) ! Calculate negative perturbation call MV_Perturb(Vars%u(i), j, -1, m%Jac%u, m%Jac%u_perturb) call HydroDyn_VarsUnpackInput(Vars, m%Jac%u_perturb, m%u_perturb) - call HydroDyn_CalcOutput(t, m%u_perturb, p, x, xd, z, OtherState, m%y_lin, m, ErrStat2, ErrMsg2); if (Failed()) return + call HydroDyn_CalcOutput(t, m%u_perturb, p, x, xd, z, OtherState, m%y_lin, m, ErrStat2, ErrMsg2, calcMorisonHstLds); if (Failed()) return call HydroDyn_VarsPackOutput(Vars, m%y_lin, m%Jac%y_neg) ! Calculate column index diff --git a/modules/hydrodyn/src/Morison.f90 b/modules/hydrodyn/src/Morison.f90 index 0b3d688142..84b51b4ff4 100644 --- a/modules/hydrodyn/src/Morison.f90 +++ b/modules/hydrodyn/src/Morison.f90 @@ -3444,7 +3444,7 @@ SUBROUTINE AllocateNodeLoadVariables(InitInp, p, m, NNodes, errStat, errMsg ) END SUBROUTINE AllocateNodeLoadVariables !---------------------------------------------------------------------------------------------------------------------------------- !> Routine for computing outputs, used in both loose and tight coupling. -SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, errMsg ) +SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, errMsg, calcHstLds ) !.................................................................................................................................. REAL(DbKi), INTENT(IN ) :: Time !< Current simulation time in seconds @@ -3459,6 +3459,8 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, TYPE(Morison_MiscVarType), INTENT(INOUT) :: m !< Misc/optimization variables INTEGER(IntKi), INTENT( OUT) :: errStat !< Error status of the operation CHARACTER(*), INTENT( OUT) :: errMsg !< Error message if errStat /= ErrID_None + LOGICAL, OPTIONAL, INTENT(IN ) :: calcHstLds !< Flag to calculate the hydrostatic loads (default: .true.) + !! Used to speed up Jacobian calculations when perturbing velocity/acceleration inputs ! Local variables INTEGER(IntKi) :: errStat2 ! Error status of the operation (occurs after initial error) @@ -3538,12 +3540,20 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, INTEGER(IntKi) :: nodeInWater REAL(SiKi) :: WaveElev1, WaveElev2, WaveElev, FDynP, FV(3), FA(3), FAMCF(3) LOGICAL :: Is1stElement, Is1stFloodedMember + LOGICAL :: calcHstLdsLocal ! Initialize errStat errStat = ErrID_None errMsg = "" Imat = 0.0_ReKi g = p%Gravity + + ! Determine whether to calculate hydrostatic loads + if (present(calcHstLds)) then + calcHstLdsLocal = calcHstLds + else + calcHstLdsLocal = .true. + end if !=============================================================================================== ! Get displaced positions of the hydrodynamic nodes @@ -3570,7 +3580,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, DO j = 1,p%FilledGroups(i)%FillNumM im = p%FilledGroups(i)%FillMList(j) IF (p%Members(im)%memfloodstatus>0) THEN - CALL getMemBallastHiPt(p%Members(im),z_hi,ErrStat2,ErrMsg2); if (Failed()) return + CALL getMemBallastHiPt(p,m,u,p%Members(im),z_hi,ErrStat2,ErrMsg2); if (Failed()) return IF ( Is1stFloodedMember ) THEN m%zFillGroup(i) = z_hi Is1stFloodedMember = .false. @@ -3658,8 +3668,6 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, ! Note: CMatrix is element local to global displaced. CTrans is the opposite. ! save some commonly used variables dl = mem%dl - z1 = pos1(3) ! get node z locations from input mesh - z2 = pos2(3) a_s1 = u%Mesh%TranslationAcc(:, mem%NodeIndx(i )) alpha_s1 = u%Mesh%RotationAcc (:, mem%NodeIndx(i )) omega_s1 = u%Mesh%RotationVel (:, mem%NodeIndx(i )) @@ -3750,17 +3758,26 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, y%Mesh%Moment(:,mem%NodeIndx(i+1)) = y%Mesh%Moment(:,mem%NodeIndx(i+1)) + F_IMG(4:6) ! ------------------- buoyancy loads: sides: Sections 3.1 and 3.2 ------------------------ - IF (mem%MHstLMod == 1) THEN + + ! Skip hydrostatic load calculation if flag is false + if (.not. calcHstLdsLocal) cycle + + ! Select hydrostatic load calculation method + select case (mem%MHstLMod) + + ! Standard hydrostatic load calculation + case (1) + IF ( p%WaveField%WaveStMod > 0_IntKi ) THEN ! If wave stretching is enabled, compute buoyancy up to free surface - CALL GetTotalWaveElev( Time, pos1, Zeta1, ErrStat2, ErrMsg2 ) - CALL GetTotalWaveElev( Time, pos2, Zeta2, ErrStat2, ErrMsg2 ) + CALL GetTotalWaveElev(p, m, Time, pos1, Zeta1, ErrStat2, ErrMsg2 ) + CALL GetTotalWaveElev(p, m, Time, pos2, Zeta2, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) ELSE ! Without wave stretching, compute buoyancy based on SWL Zeta1 = 0.0_ReKi Zeta2 = 0.0_ReKi END IF Is1stElement = ( i .EQ. 1) - CALL getElementHstLds_Mod1(mem, Time, pos1, pos2, Zeta1, Zeta2, k_hat, r1b, r2b, dl, mem%alpha(i), Is1stElement, F_B0, F_B1, F_B2, ErrStat2, ErrMsg2 ) + CALL getElementHstLds_Mod1(p, m, mem, Time, pos1, pos2, Zeta1, Zeta2, k_hat, r1b, r2b, dl, mem%alpha(i), Is1stElement, F_B0, F_B1, F_B2, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) ! Add nodal loads to mesh IF ( .NOT. Is1stElement ) THEN @@ -3774,19 +3791,24 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, y%Mesh%Moment(:,mem%NodeIndx(i )) = y%Mesh%Moment(:,mem%NodeIndx(i )) + F_B1(4:6) y%Mesh%Force (:,mem%NodeIndx(i+1)) = y%Mesh%Force (:,mem%NodeIndx(i+1)) + F_B2(1:3) y%Mesh%Moment(:,mem%NodeIndx(i+1)) = y%Mesh%Moment(:,mem%NodeIndx(i+1)) + F_B2(4:6) - ELSE IF (mem%MHstLMod == 2) THEN ! Alternative hydrostatic load calculation + + ! Alternative hydrostatic load calculation + case (2) + ! Get free surface elevation and normal at the element midpoint (both assumed constant over the element) posMid = 0.5 * (pos1+pos2) + ! rn is only used to estimate free surface normal numerically IF (mem%MSecGeom == MSecGeom_Cyl) THEN rn = 0.5 * (r1b +r2b ) ELSE IF (mem%MSecGeom == MSecGeom_Rec) THEN rn = MAX( 0.5*(Sa1b+Sa2b), 0.5*(Sb1b+Sb2b) ) END IF + IF (p%WaveField%WaveStMod > 0) THEN - CALL GetTotalWaveElev( Time, posMid, ZetaMid, ErrStat2, ErrMsg2 ) + CALL GetTotalWaveElev(p, m, Time, posMid, ZetaMid, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - CALL GetFreeSurfaceNormal( Time, posMid, rn, n_hat, ErrStat2, ErrMsg2 ) + CALL GetFreeSurfaceNormal( p, m, Time, posMid, rn, n_hat, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) FSPt = (/posMid(1),posMid(2),ZetaMid/) ! Reference point on the free surface ELSE @@ -3796,11 +3818,12 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, IF (mem%MSecGeom == MSecGeom_Cyl) THEN CALL GetSectionUnitVectors_Cyl( k_hat, y_hat, z_hat ) - CALL getElementHstLds_Mod2_Cyl( pos1, pos2, FSPt, k_hat, y_hat, z_hat, n_hat, r1b, r2b, dl, F_B1, F_B2, ErrStat2, ErrMsg2) + CALL getElementHstLds_Mod2_Cyl( p, pos1, pos2, FSPt, k_hat, y_hat, z_hat, n_hat, r1b, r2b, dl, F_B1, F_B2, ErrStat2, ErrMsg2) ELSE IF (mem%MSecGeom == MSecGeom_Rec) THEN CALL GetSectionUnitVectors_Rec( CMatrix, x_hat, y_hat ) - CALL getElementHstLds_Mod2_Rec( pos1, pos2, FSPt, k_hat, x_hat, y_hat, n_hat, Sa1b, Sa2b, Sb1b, Sb2b, dl, F_B1, F_B2, ErrStat2, ErrMsg2) + CALL getElementHstLds_Mod2_Rec( p, pos1, pos2, FSPt, k_hat, x_hat, y_hat, n_hat, Sa1b, Sa2b, Sb1b, Sb2b, dl, F_B1, F_B2, ErrStat2, ErrMsg2) END IF + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) ! Add nodal loads to mesh @@ -3810,7 +3833,8 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, y%Mesh%Moment(:,mem%NodeIndx(i )) = y%Mesh%Moment(:,mem%NodeIndx(i )) + F_B1(4:6) y%Mesh%Force (:,mem%NodeIndx(i+1)) = y%Mesh%Force (:,mem%NodeIndx(i+1)) + F_B2(1:3) y%Mesh%Moment(:,mem%NodeIndx(i+1)) = y%Mesh%Moment(:,mem%NodeIndx(i+1)) + F_B2(4:6) - END IF ! MHstLMod + + end select ! MHstLMod END DO ! i = max(mem%i_floor,1), N ! loop through member elements that are not fully buried in the seabed END IF ! NOT Modeled with Potential flow theory @@ -4039,7 +4063,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, 0.5*mem%AxCd(i)*p%WaveField%WtrDens * pi*mem%RMG(i)*dRdl_p * & ! axial part abs(dot_product( mem%k, m%vrel(:,mem%NodeIndx(i)) )) * matmul( mem%kkt, m%vrel(:,mem%NodeIndx(i)) ) ! axial part cont'd ELSE IF (mem%MSecGeom==MSecGeom_Rec) THEN - Call GetDistDrag_Rec(Time,mem,i,dSadl_p,dSbdl_p,f_hydro,ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Call GetDistDrag_Rec(p, m, u, xd, Time,mem,i,dSadl_p,dSbdl_p,f_hydro,ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) END IF CALL LumpDistrHydroLoads( f_hydro, mem%k, deltal, h_c, m%memberLoads(im)%F_D(:, i) ) y%Mesh%Force (:,mem%NodeIndx(i)) = y%Mesh%Force (:,mem%NodeIndx(i)) + m%memberLoads(im)%F_D(1:3, i) @@ -4194,7 +4218,8 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, CpFSInt = SubRatio * mem%Cp( FSElem+1) + (1.0-SubRatio) * mem%Cp( FSElem) AxCpFSInt = SubRatio * mem%AxCp(FSElem+1) + (1.0-SubRatio) * mem%AxCp(FSElem) - Call GetDistDrag_Rec(Time,mem,FSElem,dSadl_p,dSbdl_p,F_DS,ErrStat2,ErrMsg2,SubRatio,vrelFSInt); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Call GetDistDrag_Rec(p, m, u, xd, Time,mem,FSElem,dSadl_p,dSbdl_p,F_DS,ErrStat2,ErrMsg2,SubRatio,vrelFSInt) + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) ! Hydrodynamic added mass and inertia loads IF ( .NOT. mem%PropPot ) THEN @@ -4405,7 +4430,8 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, 0.5*mem%AxCd(i)*p%WaveField%WtrDens*pi*mem%RMG(i)*dRdl_p * & ! axial part abs(dot_product( mem%k, m%vrel(:,mem%NodeIndx(i)) )) * matmul( mem%kkt, m%vrel(:,mem%NodeIndx(i)) ) ! axial part cont'd ELSE IF (mem%MSecGeom==MSecGeom_Rec) THEN - Call GetDistDrag_Rec(Time,mem,i,dSadl_p,dSbdl_p,f_hydro,ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Call GetDistDrag_Rec(p, m, u, xd, Time,mem,i,dSadl_p,dSbdl_p,f_hydro,ErrStat2,ErrMsg2) + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) END IF CALL LumpDistrHydroLoads( f_hydro, mem%k, deltal, h_c, m%memberLoads(im)%F_D(:, i) ) y%Mesh%Force (:,mem%NodeIndx(i)) = y%Mesh%Force (:,mem%NodeIndx(i)) + m%memberLoads(im)%F_D(1:3, i) @@ -4574,9 +4600,9 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, if (mem%i_floor == 0) then ! both ends above or at seabed ! Compute loads on the end plate of node 1 IF (p%WaveField%WaveStMod > 0) THEN - CALL GetTotalWaveElev( Time, pos1, Zeta1, ErrStat2, ErrMsg2 ) + CALL GetTotalWaveElev(p, m, Time, pos1, Zeta1, ErrStat2, ErrMsg2) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - CALL GetFreeSurfaceNormal( Time, pos1, rn1, n_hat, ErrStat2, ErrMsg2 ) + CALL GetFreeSurfaceNormal(p, m, Time, pos1, rn1, n_hat, ErrStat2, ErrMsg2) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) FSPt = (/pos1(1),pos1(2),Zeta1/) ! Reference point on the free surface ELSE @@ -4588,7 +4614,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, CALL GetSectionUnitVectors_Cyl( k_hat1, y_hat, z_hat ) CALL GetSectionFreeSurfaceIntersects_Cyl( REAL(pos1,DbKi), REAL(FSPt,DbKi), k_hat1, y_hat, z_hat, n_hat, REAL(r1,DbKi), theta1, theta2, secStat) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - CALL GetEndPlateHstLds_Cyl(pos1, k_hat1, y_hat, z_hat, r1, theta1, theta2, F_B_End) + CALL GetEndPlateHstLds_Cyl(p, pos1, k_hat1, y_hat, z_hat, r1, theta1, theta2, F_B_End) IF (mem%MHstLMod == 1) THEN ! Check for partially wetted end plates IF ( .NOT.( EqualRealNos((theta2-theta1),0.0_DbKi) .OR. EqualRealNos((theta2-theta1),2.0_DbKi*PI_D) ) ) THEN CALL SetErrStat(ErrID_Warn, 'End plate is partially wetted with MHstLMod = 1. The buoyancy load and distribution potentially have large error. This has happened to the first node of Member ID ' //trim(num2lstr(mem%MemberID)), errStat, errMsg, RoutineName ) @@ -4596,15 +4622,15 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, END IF else if (mem%MSecGeom==MSecGeom_Rec) then CALL GetSectionUnitVectors_Rec( CMatrix1, x_hat, y_hat ) - CALL GetEndPlateHstLds_Rec(pos1, k_hat1, x_hat, y_hat, Sa1, Sb1, FSPt, n_hat, F_B_End) + CALL GetEndPlateHstLds_Rec(p, pos1, k_hat1, x_hat, y_hat, Sa1, Sb1, FSPt, n_hat, F_B_End) end if m%F_B_End(:, mem%NodeIndx( 1)) = m%F_B_End(:, mem%NodeIndx( 1)) + F_B_End ! Compute loads on the end plate of node N+1 IF (p%WaveField%WaveStMod > 0) THEN - CALL GetTotalWaveElev( Time, pos2, Zeta2, ErrStat2, ErrMsg2 ) + CALL GetTotalWaveElev(p, m, Time, pos2, Zeta2, ErrStat2, ErrMsg2) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - CALL GetFreeSurfaceNormal( Time, pos2, rn2, n_hat, ErrStat2, ErrMsg2 ) + CALL GetFreeSurfaceNormal(p, m, Time, pos2, rn2, n_hat, ErrStat2, ErrMsg2) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) FSPt = (/pos2(1),pos2(2),Zeta2/) ! Reference point on the free surface ELSE @@ -4616,7 +4642,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, CALL GetSectionUnitVectors_Cyl( k_hat2, y_hat, z_hat ) CALL GetSectionFreeSurfaceIntersects_Cyl( REAL(pos2,DbKi), REAL(FSPt,DbKi), k_hat2, y_hat, z_hat, n_hat, REAL(r2,DbKi), theta1, theta2, secStat) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - CALL GetEndPlateHstLds_Cyl(pos2, k_hat2, y_hat, z_hat, r2, theta1, theta2, F_B_End) + CALL GetEndPlateHstLds_Cyl(p, pos2, k_hat2, y_hat, z_hat, r2, theta1, theta2, F_B_End) IF (mem%MHstLMod == 1) THEN ! Check for partially wetted end plates IF ( .NOT.( EqualRealNos((theta2-theta1),0.0_DbKi) .OR. EqualRealNos((theta2-theta1),2.0_DbKi*PI_D) ) ) THEN CALL SetErrStat(ErrID_Warn, 'End plate is partially wetted with MHstLMod = 1. The buoyancy load and distribution potentially have large error. This has happened to the last node of Member ID ' //trim(num2lstr(mem%MemberID)), errStat, errMsg, RoutineName ) @@ -4624,16 +4650,16 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, END IF else if (mem%MSecGeom==MSecGeom_Rec) then CALL GetSectionUnitVectors_Rec( CMatrix2, x_hat, y_hat ) - CALL GetEndPlateHstLds_Rec(pos2, k_hat2, x_hat, y_hat, Sa2, Sb2, FSPt, n_hat, F_B_End) + CALL GetEndPlateHstLds_Rec(p, pos2, k_hat2, x_hat, y_hat, Sa2, Sb2, FSPt, n_hat, F_B_End) end if m%F_B_End(:, mem%NodeIndx(N+1)) = m%F_B_End(:, mem%NodeIndx(N+1)) - F_B_End elseif ( mem%doEndBuoyancy ) then ! The member crosses the seabed line so only the upper end potentially have hydrostatic load ! Only compute the loads on the end plate of node N+1 IF (p%WaveField%WaveStMod > 0) THEN - CALL GetTotalWaveElev( Time, pos2, Zeta2, ErrStat2, ErrMsg2 ) + CALL GetTotalWaveElev(p, m, Time, pos2, Zeta2, ErrStat2, ErrMsg2) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - CALL GetFreeSurfaceNormal( Time, pos2, rn2, n_hat, ErrStat2, ErrMsg2 ) + CALL GetFreeSurfaceNormal(p, m, Time, pos2, rn2, n_hat, ErrStat2, ErrMsg2) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) FSPt = (/pos2(1),pos2(2),Zeta2/) ! Reference point on the free surface ELSE @@ -4645,7 +4671,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, CALL GetSectionUnitVectors_Cyl( k_hat2, y_hat, z_hat ) CALL GetSectionFreeSurfaceIntersects_Cyl( REAL(pos2,DbKi), REAL(FSPt,DbKi), k_hat2, y_hat, z_hat, n_hat, REAL(r2,DbKi), theta1, theta2, secStat) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - CALL GetEndPlateHstLds_Cyl(pos2, k_hat2, y_hat, z_hat, r2, theta1, theta2, F_B_End) + CALL GetEndPlateHstLds_Cyl(p, pos2, k_hat2, y_hat, z_hat, r2, theta1, theta2, F_B_End) IF (mem%MHstLMod == 1) THEN ! Check for partially wetted end plates IF ( .NOT.( EqualRealNos((theta2-theta1),0.0_DbKi) .OR. EqualRealNos((theta2-theta1),2.0_DbKi*PI_D) ) ) THEN CALL SetErrStat(ErrID_Warn, 'End plate is partially wetted with MHstLMod = 1. The buoyancy load and distribution potentially have large error. This has happened to the last node of Member ID ' //trim(num2lstr(mem%MemberID)), errStat, errMsg, RoutineName ) @@ -4653,7 +4679,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, END IF else if (mem%MSecGeom==MSecGeom_Rec) then CALL GetSectionUnitVectors_Rec( CMatrix2, x_hat, y_hat ) - CALL GetEndPlateHstLds_Rec(pos2, k_hat2, x_hat, y_hat, Sa2, Sb2, FSPt, n_hat, F_B_End) + CALL GetEndPlateHstLds_Rec(p, pos2, k_hat2, x_hat, y_hat, Sa2, Sb2, FSPt, n_hat, F_B_End) end if m%F_B_End(:, mem%NodeIndx(N+1)) = m%F_B_End(:, mem%NodeIndx(N+1)) - F_B_End @@ -4687,7 +4713,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, ! Effect of wave stretching already baked into m%FDynP, m%FA, and m%vrel. No additional modification needed. ! Joint yaw offset - call YawJoint(J,u%PtfmRefY,AM_End,An_End,DP_Const_End,I_MG_End,ErrStat2,ErrMsg2) + call YawJoint(p, J,u%PtfmRefY,AM_End,An_End,DP_Const_End,I_MG_End,ErrStat2,ErrMsg2) call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) ! Lumped added mass loads @@ -4750,7 +4776,19 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, CONTAINS - SUBROUTINE GetTotalWaveElev( Time, pos, Zeta, ErrStat, ErrMsg ) + logical function Failed() + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Failed = ErrStat >= AbortErrLev + !if (Failed) then + ! call FailCleanup() + !endif + end function Failed + +END SUBROUTINE Morison_CalcOutput + + SUBROUTINE GetTotalWaveElev(p, m, Time, pos, Zeta, ErrStat, ErrMsg ) + TYPE(Morison_ParameterType), INTENT( IN ) :: p + TYPE(Morison_MiscVarType), INTENT( INOUT ) :: m REAL(DbKi), INTENT( IN ) :: Time REAL(ReKi), INTENT( IN ) :: pos(*) ! Position at which free-surface elevation is to be calculated. Third entry ignored if present. REAL(ReKi), INTENT( OUT ) :: Zeta ! Total free-surface elevation with first- and second-order contribution (if present) @@ -4767,7 +4805,9 @@ SUBROUTINE GetTotalWaveElev( Time, pos, Zeta, ErrStat, ErrMsg ) END SUBROUTINE GetTotalWaveElev - SUBROUTINE GetFreeSurfaceNormal( Time, pos, r, n, ErrStat, ErrMsg) + SUBROUTINE GetFreeSurfaceNormal(p, m, Time, pos, r, n, ErrStat, ErrMsg) + TYPE(Morison_ParameterType), INTENT( IN ) :: p + TYPE(Morison_MiscVarType), INTENT( INOUT ) :: m REAL(DbKi), INTENT( In ) :: Time REAL(ReKi), INTENT( In ) :: pos(:) ! Position at which free-surface normal is to be calculated. Third entry ignored if present. REAL(ReKi), INTENT( In ) :: r ! Distance for central differencing @@ -4859,8 +4899,8 @@ SUBROUTINE GetSectionFreeSurfaceIntersects_Cyl( pos0, FSPt, k_hat, y_hat, z_hat, END SUBROUTINE GetSectionFreeSurfaceIntersects_Cyl - SUBROUTINE GetSectionHstLds_Cyl( origin, pos0, k_hat, y_hat, z_hat, R, dRdl, theta1, theta2, dFdl) - + SUBROUTINE GetSectionHstLds_Cyl(p, origin, pos0, k_hat, y_hat, z_hat, R, dRdl, theta1, theta2, dFdl) + TYPE(Morison_ParameterType), INTENT( IN ) :: p REAL(DbKi), INTENT( IN ) :: origin(3) REAL(DbKi), INTENT( IN ) :: pos0(3) REAL(DbKi), INTENT( IN ) :: k_hat(3) @@ -4888,12 +4928,12 @@ SUBROUTINE GetSectionHstLds_Cyl( origin, pos0, k_hat, y_hat, z_hat, R, dRdl, the dFdl(1:3) = -R *dRdl*C0*k_hat + R*C1*y_hat + R*C2*z_hat dFdl(4:6) = -R**2*dRdl*C2*y_hat + R**2*dRdl*C1*z_hat + CROSS_PRODUCT((pos0-origin),dFdl(1:3)) - dFdl = dFdl * p%WaveField%WtrDens * g + dFdl = dFdl * p%WaveField%WtrDens * p%Gravity END SUBROUTINE GetSectionHstLds_Cyl - SUBROUTINE GetSectionHstLds_Rec( origin, pos0, k_hat, x_hat, y_hat, Sa, Sb, dSadl, dSbdl, rFS, nFS, dFdl, secStat) - + SUBROUTINE GetSectionHstLds_Rec(p, origin, pos0, k_hat, x_hat, y_hat, Sa, Sb, dSadl, dSbdl, rFS, nFS, dFdl, secStat) + TYPE(Morison_ParameterType), INTENT( IN ) :: p REAL(DbKi), INTENT( IN ) :: origin(3) REAL(DbKi), INTENT( IN ) :: pos0(3) REAL(DbKi), INTENT( IN ) :: k_hat(3) @@ -4996,12 +5036,12 @@ SUBROUTINE GetSectionHstLds_Rec( origin, pos0, k_hat, x_hat, y_hat, Sa, Sb, dSad end do - dFdl = dFdl * p%WaveField%WtrDens * g + dFdl = dFdl * p%WaveField%WtrDens * p%Gravity END SUBROUTINE GetSectionHstLds_Rec - SUBROUTINE getElementHstLds_Mod2_Cyl( pos1In, pos2In, FSPtIn, k_hatIn, y_hatIn, z_hatIn, n_hatIn, r1In, r2In, dlIn, F_B1, F_B2, ErrStat, ErrMsg ) - + SUBROUTINE getElementHstLds_Mod2_Cyl(p, pos1In, pos2In, FSPtIn, k_hatIn, y_hatIn, z_hatIn, n_hatIn, r1In, r2In, dlIn, F_B1, F_B2, ErrStat, ErrMsg ) + TYPE(Morison_ParameterType), INTENT( IN ) :: p REAL(ReKi), INTENT( IN ) :: pos1In(3) REAL(ReKi), INTENT( IN ) :: pos2In(3) REAL(ReKi), INTENT( IN ) :: FSPtIn(3) @@ -5024,7 +5064,6 @@ SUBROUTINE getElementHstLds_Mod2_Cyl( pos1In, pos2In, FSPtIn, k_hatIn, y_hatIn, INTEGER(IntKi) :: errStat2 CHARACTER(ErrMsgLen) :: errMsg2 ErrStat = ErrID_None - ErrMsg = "" pos1 = REAL(pos1In,DbKi) pos2 = REAL(pos2In,DbKi) @@ -5055,19 +5094,21 @@ SUBROUTINE getElementHstLds_Mod2_Cyl( pos1In, pos2In, FSPtIn, k_hatIn, y_hatIn, ! Section load at node 1 CALL GetSectionFreeSurfaceIntersects_Cyl( pos1, FSPt, REAL(k_hat,ReKi), REAL(y_hat,ReKi), REAL(z_hat,ReKi), REAL(n_hat,ReKi), r1, theta1, theta2, secStat1) - CALL GetSectionHstLds_Cyl( pos1, pos1, k_hat, y_hat, z_hat, r1, dRdl, theta1, theta2, dFdl1) + CALL GetSectionHstLds_Cyl(p, pos1, pos1, k_hat, y_hat, z_hat, r1, dRdl, theta1, theta2, dFdl1) ! Section load at midpoint CALL GetSectionFreeSurfaceIntersects_Cyl( posMid, FSPt, REAL(k_hat,ReKi), REAL(y_hat,ReKi), REAL(z_hat,ReKi), REAL(n_hat,ReKi), rMid, theta1, theta2, secStatMid) - CALL GetSectionHstLds_Cyl( pos1, posMid, k_hat, y_hat, z_hat, rMid, dRdl, theta1, theta2, dFdlMid) + CALL GetSectionHstLds_Cyl(p, pos1, posMid, k_hat, y_hat, z_hat, rMid, dRdl, theta1, theta2, dFdlMid) ! Section load at node 2 CALL GetSectionFreeSurfaceIntersects_Cyl( pos2, FSPt, REAL(k_hat,ReKi), REAL(y_hat,ReKi), REAL(z_hat,ReKi), REAL(n_hat,ReKi), r2, theta1, theta2, secStat2) - CALL GetSectionHstLds_Cyl( pos1, pos2, k_hat, y_hat, z_hat, r2, dRdl, theta1, theta2, dFdl2) + CALL GetSectionHstLds_Cyl(p, pos1, pos2, k_hat, y_hat, z_hat, r2, dRdl, theta1, theta2, dFdl2) ! Adaptively refine the load integration over the element - CALL RefineElementHstLds_Cyl(pos1,pos1,posMid,pos2,FSPt,r1,rMid,r2,dl,dRdl,secStat1,secStatMid,secStat2,k_hat,y_hat,z_hat,n_hat,dFdl1,dFdlMid,dFdl2,1,F_B,ErrStat2,ErrMsg2) - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + CALL RefineElementHstLds_Cyl(p,pos1,pos1,posMid,pos2,FSPt,r1,rMid,r2,dl,dRdl,secStat1,secStatMid,secStat2,k_hat,y_hat,z_hat,n_hat,dFdl1,dFdlMid,dFdl2,1,F_B,ErrStat2) + if (ErrStat2 /= ErrID_None) then + CALL SetErrStat( ErrStat2, 'Tolerance for element hydrostatic load not met after the maximum allowed level of recursion is reached. Consider reducing MDivSize.', ErrStat, ErrMsg, RoutineName ) + end if ! Distribute the hydrostatic load to the two end nodes F_B1(1:3) = 0.5_DbKi * F_B(1:3) @@ -5078,8 +5119,8 @@ SUBROUTINE getElementHstLds_Mod2_Cyl( pos1In, pos2In, FSPtIn, k_hatIn, y_hatIn, END SUBROUTINE getElementHstLds_Mod2_Cyl - SUBROUTINE getElementHstLds_Mod2_Rec( pos1In, pos2In, FSPtIn, k_hatIn, x_hatIn, y_hatIn, n_hatIn, Sa1In, Sa2In, Sb1In, Sb2In, dlIn, F_B1, F_B2, ErrStat, ErrMsg ) - + SUBROUTINE getElementHstLds_Mod2_Rec(p, pos1In, pos2In, FSPtIn, k_hatIn, x_hatIn, y_hatIn, n_hatIn, Sa1In, Sa2In, Sb1In, Sb2In, dlIn, F_B1, F_B2, ErrStat, ErrMsg ) + TYPE(Morison_ParameterType), INTENT( IN ) :: p REAL(ReKi), INTENT( IN ) :: pos1In(3) REAL(ReKi), INTENT( IN ) :: pos2In(3) REAL(ReKi), INTENT( IN ) :: FSPtIn(3) @@ -5136,16 +5177,16 @@ SUBROUTINE getElementHstLds_Mod2_Rec( pos1In, pos2In, FSPtIn, k_hatIn, x_hatIn, END IF ! Section load at node 1 - CALL GetSectionHstLds_Rec( pos1, pos1, k_hat, x_hat, y_hat, Sa1, Sb1, dSadl, dSbdl, FSPt, n_hat, dFdl1, secStat1) + CALL GetSectionHstLds_Rec(p, pos1, pos1, k_hat, x_hat, y_hat, Sa1, Sb1, dSadl, dSbdl, FSPt, n_hat, dFdl1, secStat1) ! Section load at midpoint - CALL GetSectionHstLds_Rec( pos1, posMid, k_hat, x_hat, y_hat, SaMid, SbMid, dSadl, dSbdl, FSPt, n_hat, dFdlMid, secStatMid) + CALL GetSectionHstLds_Rec(p, pos1, posMid, k_hat, x_hat, y_hat, SaMid, SbMid, dSadl, dSbdl, FSPt, n_hat, dFdlMid, secStatMid) ! Section load at node 2 - CALL GetSectionHstLds_Rec( pos1, pos2, k_hat, x_hat, y_hat, Sa2, Sb2, dSadl, dSbdl, FSPt, n_hat, dFdl2, secStat2) + CALL GetSectionHstLds_Rec(p, pos1, pos2, k_hat, x_hat, y_hat, Sa2, Sb2, dSadl, dSbdl, FSPt, n_hat, dFdl2, secStat2) ! Adaptively refine the load integration over the element - CALL RefineElementHstLds_Rec(pos1,pos1,posMid,pos2,FSPt,Sa1,SaMid,Sa2,Sb1,SbMid,Sb2,dl,dSadl,dSbdl, & + CALL RefineElementHstLds_Rec(p,pos1,pos1,posMid,pos2,FSPt,Sa1,SaMid,Sa2,Sb1,SbMid,Sb2,dl,dSadl,dSbdl, & secStat1,secStatMid,secStat2,k_hat,x_hat,y_hat,n_hat,dFdl1,dFdlMid,dFdl2,1,F_B,ErrStat2,ErrMsg2) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) @@ -5158,8 +5199,8 @@ SUBROUTINE getElementHstLds_Mod2_Rec( pos1In, pos2In, FSPtIn, k_hatIn, x_hatIn, END SUBROUTINE getElementHstLds_Mod2_Rec - RECURSIVE SUBROUTINE RefineElementHstLds_Cyl( origin, pos1, posMid, pos2, FSPt, r1, rMid, r2, dl, dRdl,secStat1,secStatMid,secStat2, k_hat, y_hat, z_hat, n_hat, dFdl1, dFdlMid, dFdl2, recurLvl, F_B_5pt, ErrStat, ErrMsg) - + RECURSIVE SUBROUTINE RefineElementHstLds_Cyl(p, origin, pos1, posMid, pos2, FSPt, r1, rMid, r2, dl, dRdl,secStat1,secStatMid,secStat2, k_hat, y_hat, z_hat, n_hat, dFdl1, dFdlMid, dFdl2, recurLvl, F_B_5pt, ErrStat) + TYPE(Morison_ParameterType), INTENT( IN ) :: p REAL(DbKi), INTENT( IN ) :: origin(3) REAL(DbKi), INTENT( IN ) :: pos1(3) REAL(DbKi), INTENT( IN ) :: posMid(3) @@ -5182,15 +5223,15 @@ RECURSIVE SUBROUTINE RefineElementHstLds_Cyl( origin, pos1, posMid, pos2, FSPt, REAL(DbKi), INTENT( IN ) :: dFdl2(6) INTEGER(IntKi), INTENT( IN ) :: recurLvl REAL(DbKi), INTENT( OUT ) :: F_B_5pt(6) + INTEGER(IntKi), INTENT( OUT ) :: ErrStat ! Error status of the operation - CHARACTER(*), INTENT( OUT ) :: ErrMsg ! Error message if errStat /= ErrID_None REAL(DbKi) :: theta1,theta2 REAL(DbKi) :: posMidL(3), posMidR(3) REAL(DbKi) :: rMidL, rMidR REAL(DbKi) :: dFdlMidL(6), dFdlMidR(6), F_B_3pt(6) REAL(DbKi) :: error(6), tmp(6) - LOGICAL :: refine, tolMet + LOGICAL :: tolMet INTEGER(IntKi) :: i INTEGER(IntKi) :: secStatMidL, secStatMidR REAL(DbKi), PARAMETER :: RelTol = 1.0E-6 @@ -5199,7 +5240,6 @@ RECURSIVE SUBROUTINE RefineElementHstLds_Cyl( origin, pos1, posMid, pos2, FSPt, CHARACTER(*), PARAMETER :: RoutineName = "RefineElementHstLds_Cyl" ErrStat = ErrID_None - ErrMsg = "" posMidL = 0.5_DbKi*(pos1+posMid) posMidR = 0.5_DbKi*(posMid+pos2) @@ -5221,45 +5261,37 @@ RECURSIVE SUBROUTINE RefineElementHstLds_Cyl( origin, pos1, posMid, pos2, FSPt, ! Mid point of left section CALL GetSectionFreeSurfaceIntersects_Cyl( posMidL, FSPt, REAL(k_hat,ReKi), REAL(y_hat,ReKi), REAL(z_hat,ReKi), REAL(n_hat,ReKi), rMidL, theta1, theta2, secStatMidL) - CALL GetSectionHstLds_Cyl( origin, posMidL, k_hat, y_hat, z_hat, rMidL, dRdl, theta1, theta2, dFdlMidL) + CALL GetSectionHstLds_Cyl(p, origin, posMidL, k_hat, y_hat, z_hat, rMidL, dRdl, theta1, theta2, dFdlMidL) ! Mid point of right section CALL GetSectionFreeSurfaceIntersects_Cyl( posMidR, FSPt, REAL(k_hat,ReKi), REAL(y_hat,ReKi), REAL(z_hat,ReKi), REAL(n_hat,ReKi), rMidR, theta1, theta2, secStatMidR) - CALL GetSectionHstLds_Cyl( origin, posMidR, k_hat, y_hat, z_hat, rMidR, dRdl, theta1, theta2, dFdlMidR) + CALL GetSectionHstLds_Cyl(p, origin, posMidR, k_hat, y_hat, z_hat, rMidR, dRdl, theta1, theta2, dFdlMidR) F_B_5pt = (dFdl1 + 4.0_DbKi*dFdlMidL + 2.0_DbKi*dFdlMid + 4.0_DbKi*dFdlMidR + dFdl2) * dl/12.0_DbKi + ! Calculate error and check against tolerance error = ABS(F_B_3pt - F_B_5pt) - tolMet = .TRUE. - DO i = 1,6 - IF ( error(i) > MAX(RelTol*ABS(F_B_5pt(i)),AbsTol) ) THEN - tolMet = .FALSE. - END IF - END DO - refine = .NOT. tolMet - IF (ABS(secStat1-secStat2)>1) THEN ! (Sub)element bounds the waterplane - refine = .TRUE. ! Keep refining irrespective of tolMet to avoid premature termination - END IF - IF ( recurLvl > maxRecurLvl ) THEN - refine = .FALSE. - IF (.NOT. tolMet) THEN - CALL SetErrStat(ErrID_Warn, 'Tolerance for element hydrostatic load not met after the maximum allowed level of recursion is reached. Consider reducing MDivSize.', ErrStat, ErrMsg, RoutineName ) - ! ELSE - ! Free surface is likely normal to the element. - END IF - END IF - - IF (refine) THEN ! Recursively refine the load integration if tolerance not met - CALL RefineElementHstLds_Cyl(origin,pos1,posMidL,posMid,FSPt,r1,rMidL,rMid,0.5_DbKi*dl,dRdl,secStat1,secStatMidL,secStatMid,k_hat,y_hat,z_hat,n_hat,dFdl1,dFdlMidL,dFdlMid, recurLvl+1, tmp, ErrStat, ErrMsg) - CALL RefineElementHstLds_Cyl(origin,posMid,posMidR,pos2,FSPt,rMid,rMidR,r2,0.5_DbKi*dl,dRdl,secStatMid,secStatMidR,secStat2,k_hat,y_hat,z_hat,n_hat,dFdlMid,dFdlMidR,dFdl2, recurLvl+1, F_B_5pt, ErrStat, ErrMsg) - F_B_5pt = F_B_5pt + tmp - END IF + tolMet = all(error <= MAX(RelTol*ABS(F_B_5pt),AbsTol)) + + ! If tolerance was met and (sub)element does not bound the waterplane, exit recursion + if (tolMet .and. (ABS(secStat1-secStat2) <= 1)) return + + ! If maximum recursion level reached, exit recursion with warning + if (recurLvl > maxRecurLvl) then + ErrStat = ErrID_Warn + return + end if + ! Recursively refine the load integration if tolerance not met + CALL RefineElementHstLds_Cyl(p,origin,pos1,posMidL,posMid,FSPt,r1,rMidL,rMid,0.5_DbKi*dl,dRdl,secStat1,secStatMidL,secStatMid,k_hat,y_hat,z_hat,n_hat,dFdl1,dFdlMidL,dFdlMid, recurLvl+1, tmp, ErrStat) + CALL RefineElementHstLds_Cyl(p,origin,posMid,posMidR,pos2,FSPt,rMid,rMidR,r2,0.5_DbKi*dl,dRdl,secStatMid,secStatMidR,secStat2,k_hat,y_hat,z_hat,n_hat,dFdlMid,dFdlMidR,dFdl2, recurLvl+1, F_B_5pt, ErrStat) + F_B_5pt = F_B_5pt + tmp + END SUBROUTINE RefineElementHstLds_Cyl - RECURSIVE SUBROUTINE RefineElementHstLds_Rec( origin, pos1, posMid, pos2, FSPt, Sa1, SaMid, Sa2, Sb1, SbMid, Sb2, dl, dSadl, dSbdl, & + RECURSIVE SUBROUTINE RefineElementHstLds_Rec(p, origin, pos1, posMid, pos2, FSPt, Sa1, SaMid, Sa2, Sb1, SbMid, Sb2, dl, dSadl, dSbdl, & secStat1, secStatMid, secStat2, k_hat, x_hat, y_hat, n_hat, dFdl1, dFdlMid, dFdl2, recurLvl, F_B_5pt, ErrStat, ErrMsg) - + TYPE(Morison_ParameterType), INTENT( IN ) :: p REAL(DbKi), INTENT( IN ) :: origin(3) REAL(DbKi), INTENT( IN ) :: pos1(3) REAL(DbKi), INTENT( IN ) :: posMid(3) @@ -5321,10 +5353,10 @@ RECURSIVE SUBROUTINE RefineElementHstLds_Rec( origin, pos1, posMid, pos2, FSPt, F_B_3pt = (dFdl1 + 4.0_DbKi*dFdlMid + dFdl2) * dl/6.0_DbKi ! Mid point of left section - CALL GetSectionHstLds_Rec( origin, posMidL, k_hat, x_hat, y_hat, SaMidL, SbMidL, dSadl, dSbdl, FSPt, n_hat, dFdlMidL, secStatMidL) + CALL GetSectionHstLds_Rec(p, origin, posMidL, k_hat, x_hat, y_hat, SaMidL, SbMidL, dSadl, dSbdl, FSPt, n_hat, dFdlMidL, secStatMidL) ! Mid point of right section - CALL GetSectionHstLds_Rec( origin, posMidR, k_hat, x_hat, y_hat, SaMidR, SbMidR, dSadl, dSbdl, FSPt, n_hat, dFdlMidR, secStatMidR) + CALL GetSectionHstLds_Rec(p, origin, posMidR, k_hat, x_hat, y_hat, SaMidR, SbMidR, dSadl, dSbdl, FSPt, n_hat, dFdlMidR, secStatMidR) F_B_5pt = (dFdl1 + 4.0_DbKi*dFdlMidL + 2.0_DbKi*dFdlMid + 4.0_DbKi*dFdlMidR + dFdl2) * dl/12.0_DbKi @@ -5349,17 +5381,17 @@ RECURSIVE SUBROUTINE RefineElementHstLds_Rec( origin, pos1, posMid, pos2, FSPt, END IF IF (refine) THEN ! Recursively refine the load integration if tolerance not met - CALL RefineElementHstLds_Rec(origin,pos1,posMidL,posMid,FSPt,Sa1,SaMidL,SaMid,Sb1,SbMidL,SbMid,0.5_DbKi*dl,dSadl,dSbdl, & + CALL RefineElementHstLds_Rec(p,origin,pos1,posMidL,posMid,FSPt,Sa1,SaMidL,SaMid,Sb1,SbMidL,SbMid,0.5_DbKi*dl,dSadl,dSbdl, & secStat1,secStatMidL,secStatMid,k_hat,x_hat,y_hat,n_hat,dFdl1,dFdlMidL,dFdlMid,recurLvl+1,tmp,ErrStat,ErrMsg) - CALL RefineElementHstLds_Rec(origin,posMid,posMidR,pos2,FSPt,SaMid,SaMidR,Sa2,SbMid,SbMidR,Sb2,0.5_DbKi*dl,dSadl,dSbdl, & + CALL RefineElementHstLds_Rec(p,origin,posMid,posMidR,pos2,FSPt,SaMid,SaMidR,Sa2,SbMid,SbMidR,Sb2,0.5_DbKi*dl,dSadl,dSbdl, & secStatMid,secStatMidR,secStat2,k_hat,x_hat,y_hat,n_hat,dFdlMid,dFdlMidR,dFdl2,recurLvl+1,F_B_5pt,ErrStat,ErrMsg) F_B_5pt = F_B_5pt + tmp END IF END SUBROUTINE RefineElementHstLds_Rec - SUBROUTINE GetEndPlateHstLds_Cyl(pos0, k_hat, y_hat, z_hat, R, theta1, theta2, F) - + SUBROUTINE GetEndPlateHstLds_Cyl(p, pos0, k_hat, y_hat, z_hat, R, theta1, theta2, F) + TYPE(Morison_ParameterType), INTENT( IN ) :: p REAL(ReKi), INTENT( IN ) :: pos0(3) REAL(ReKi), INTENT( IN ) :: k_hat(3) REAL(ReKi), INTENT( IN ) :: y_hat(3) @@ -5403,7 +5435,7 @@ SUBROUTINE GetEndPlateHstLds_Cyl(pos0, k_hat, y_hat, z_hat, R, theta1, theta2, F ! End plate force in the k_hat direction Fk = -0.5_DbKi*Z0*(R_2*dTheta-tmp1) + cosPhi/6.0_DbKi*( 2.0_DbKi*dy_3 - z1*z2*dy - z1_2*(y2+2.0_DbKi*y1) + z2_2*(y1+2.0_DbKi*y2) ) - F(1:3) = p%WaveField%WtrDens * g * Fk * k_hat + F(1:3) = p%WaveField%WtrDens * p%Gravity * Fk * k_hat ! End plate moment in the y_hat and z_hat direction My = Z0/6.0_DbKi*( 2.0_DbKi*dy_3 + 2.0_DbKi*dy*tmp2 + 3.0_DbKi*tmp1*sz ) & ! y_hat component @@ -5425,12 +5457,13 @@ SUBROUTINE GetEndPlateHstLds_Cyl(pos0, k_hat, y_hat, z_hat, R, theta1, theta2, F Mz = -Z0/ 6.0_DbKi*( dz*( y2*y2 + y2*y1 + y1*y1 + tmp2 - 3.0_DbKi*R_2 ) ) & -cosPhi/24.0_DbKi*( dz_2*(3.0_DbKi*z2_2+3.0_DbKi*z1_2+2.0_DbKi*y1*y2-6.0_DbKi*R_2) + dz*dz*(y1*y1-y2*y2) + 4.0_DbKi*dz*(y1*y1*z1+y2*y2*z2) ) - F(4:6) = p%WaveField%WtrDens * g * (My*y_hat + Mz*z_hat) + F(4:6) = p%WaveField%WtrDens * p%Gravity * (My*y_hat + Mz*z_hat) END SUBROUTINE GetEndPlateHstLds_Cyl - SUBROUTINE GetEndPlateHstLds_Rec(pos0, k_hat, x_hat, y_hat, Sa, Sb, rFS, nFS, F) + SUBROUTINE GetEndPlateHstLds_Rec(p, pos0, k_hat, x_hat, y_hat, Sa, Sb, rFS, nFS, F) + TYPE(Morison_ParameterType), INTENT( IN ) :: p REAL(ReKi), INTENT( IN ) :: pos0(3) REAL(ReKi), INTENT( IN ) :: k_hat(3) REAL(ReKi), INTENT( IN ) :: x_hat(3) @@ -5497,7 +5530,7 @@ SUBROUTINE GetEndPlateHstLds_Rec(pos0, k_hat, x_hat, y_hat, Sa, Sb, rFS, nFS, F) s1 = -0.5*Sa s2 = Sa * (0.5 - dot_product(rFS-rv(:,3),nFS)/dot_product(rv(:,4)-rv(:,3),nFS) ) end if - call GetHstLdsOnTrapezoid(REAL(pos0,DbKi),s1,s2,h1s1,h1s2,h2s1,h2s2,REAL(k_hat,DbKi),REAL(x_hat,DbKi),REAL(y_hat,DbKi),Ftmp1) + call GetHstLdsOnTrapezoid(p, REAL(pos0,DbKi),s1,s2,h1s1,h1s2,h2s1,h2s2,REAL(k_hat,DbKi),REAL(x_hat,DbKi),REAL(y_hat,DbKi),Ftmp1) F = Ftmp1 else if (numVInWtr == 2) then ! Two neighboring vertices in water if (vInWtr(1) .and. vInWtr(2)) then @@ -5520,7 +5553,7 @@ SUBROUTINE GetEndPlateHstLds_Rec(pos0, k_hat, x_hat, y_hat, Sa, Sb, rFS, nFS, F) h1s2 = -0.5*Sb s1 = Sa * (-0.5 + dot_product(rFS-rv(:,1),nFS)/dot_product(rv(:,2)-rv(:,1),nFS) ) s2 = Sa * (0.5 - dot_product(rFS-rv(:,3),nFS)/dot_product(rv(:,4)-rv(:,3),nFS) ) - call GetHstLdsOnTrapezoid(REAL(pos0,DbKi),s2,0.5_DbKi*Sa,-0.5_DbKi*Sb,-0.5_DbKi*Sb,0.5_DbKi*Sb,0.5_DbKi*Sb,REAL(k_hat,DbKi),REAL(x_hat,DbKi),REAL(y_hat,DbKi),Ftmp2) + call GetHstLdsOnTrapezoid(p, REAL(pos0,DbKi),s2,0.5_DbKi*Sa,-0.5_DbKi*Sb,-0.5_DbKi*Sb,0.5_DbKi*Sb,0.5_DbKi*Sb,REAL(k_hat,DbKi),REAL(x_hat,DbKi),REAL(y_hat,DbKi),Ftmp2) else if (vInWtr(3) .and. vInWtr(4)) then ! Sides 2 & 4 intersects the free surface ! Side 3 submerged and side 1 dry @@ -5541,9 +5574,9 @@ SUBROUTINE GetEndPlateHstLds_Rec(pos0, k_hat, x_hat, y_hat, Sa, Sb, rFS, nFS, F) h1s2 = -0.5*Sb s1 = Sa * ( 0.5 - dot_product(rFS-rv(:,3),nFS)/dot_product(rv(:,4)-rv(:,3),nFS) ) s2 = Sa * (-0.5 + dot_product(rFS-rv(:,1),nFS)/dot_product(rv(:,2)-rv(:,1),nFS) ) - call GetHstLdsOnTrapezoid(REAL(pos0,DbKi),-0.5_DbKi*Sa,s1,-0.5_DbKi*Sb,-0.5_DbKi*Sb,0.5_DbKi*Sb,0.5_DbKi*Sb,REAL(k_hat,DbKi),REAL(x_hat,DbKi),REAL(y_hat,DbKi),Ftmp2) + call GetHstLdsOnTrapezoid(p,REAL(pos0,DbKi),-0.5_DbKi*Sa,s1,-0.5_DbKi*Sb,-0.5_DbKi*Sb,0.5_DbKi*Sb,0.5_DbKi*Sb,REAL(k_hat,DbKi),REAL(x_hat,DbKi),REAL(y_hat,DbKi),Ftmp2) end if - call GetHstLdsOnTrapezoid(REAL(pos0,DbKi),s1,s2,h1s1,h1s2,h2s1,h2s2,REAL(k_hat,DbKi),REAL(x_hat,DbKi),REAL(y_hat,DbKi),Ftmp1) + call GetHstLdsOnTrapezoid(p,REAL(pos0,DbKi),s1,s2,h1s1,h1s2,h2s1,h2s2,REAL(k_hat,DbKi),REAL(x_hat,DbKi),REAL(y_hat,DbKi),Ftmp1) F = Ftmp1 + Ftmp2 else if (numVInWtr == 3) then ! Only one vertex out of water if (.not. vInWtr(1)) then @@ -5579,17 +5612,18 @@ SUBROUTINE GetEndPlateHstLds_Rec(pos0, k_hat, x_hat, y_hat, Sa, Sb, rFS, nFS, F) s1 = -0.5*Sa s2 = Sa * ( 0.5 - dot_product(rFS-rv(:,3),nFS)/dot_product(rv(:,4)-rv(:,3),nFS) ) end if - call GetHstLdsOnTrapezoid(REAL(pos0,DbKi),s1,s2,h1s1,h1s2,h2s1,h2s2,REAL(k_hat,DbKi),REAL(x_hat,DbKi),REAL(y_hat,DbKi),Ftmp1) - F(1:3) = -p%WaveField%WtrDens*g*z0*Sa*Sb*k_hat - Ftmp1(1:3) - F(4:6) = p%WaveField%WtrDens*g*(Sa**3*Sb*x_hat(3)*y_hat-Sa*Sb**3*y_hat(3)*x_hat)/12.0 - Ftmp1(4:6) + call GetHstLdsOnTrapezoid(p,REAL(pos0,DbKi),s1,s2,h1s1,h1s2,h2s1,h2s2,REAL(k_hat,DbKi),REAL(x_hat,DbKi),REAL(y_hat,DbKi),Ftmp1) + F(1:3) = -p%WaveField%WtrDens*p%Gravity*z0*Sa*Sb*k_hat - Ftmp1(1:3) + F(4:6) = p%WaveField%WtrDens*p%Gravity*(Sa**3*Sb*x_hat(3)*y_hat-Sa*Sb**3*y_hat(3)*x_hat)/12.0 - Ftmp1(4:6) else if (numVInWtr == 4) then ! Submerged endplate - F(1:3) = -p%WaveField%WtrDens*g*z0*Sa*Sb*k_hat - F(4:6) = p%WaveField%WtrDens*g*(Sa**3*Sb*x_hat(3)*y_hat-Sa*Sb**3*y_hat(3)*x_hat)/12.0 + F(1:3) = -p%WaveField%WtrDens*p%Gravity*z0*Sa*Sb*k_hat + F(4:6) = p%WaveField%WtrDens*p%Gravity*(Sa**3*Sb*x_hat(3)*y_hat-Sa*Sb**3*y_hat(3)*x_hat)/12.0 end if END SUBROUTINE GetEndPlateHstLds_Rec - SUBROUTINE GetHstLdsOnTrapezoid(pos0, s1, s2, h1s1, h1s2, h2s1, h2s2, k_hat, x_hat, y_hat, F) + SUBROUTINE GetHstLdsOnTrapezoid(p, pos0, s1, s2, h1s1, h1s2, h2s1, h2s2, k_hat, x_hat, y_hat, F) + TYPE(Morison_ParameterType), INTENT( IN ) :: p REAL(DbKi), INTENT( IN ) :: pos0(3) REAL(DbKi), INTENT( IN ) :: s1, s2, h1s1, h1s2, h2s1, h2s2 REAL(DbKi), INTENT( IN ) :: k_hat(3), x_hat(3), y_hat(3) @@ -5639,12 +5673,14 @@ SUBROUTINE GetHstLdsOnTrapezoid(pos0, s1, s2, h1s1, h1s2, h2s1, h2s2, k_hat, x_h +x_hat(3)/12.0*(3.0*dp*ds4+4.0*dq*ds3) & +y_hat(3)/24.0*tmp )*y_hat - F = p%WaveField%WtrDens * g * F + F = p%WaveField%WtrDens * p%Gravity * F END SUBROUTINE GetHstLdsOnTrapezoid - SUBROUTINE getElementHstLds_Mod1( mem, Time, pos1, pos2, Zeta1, Zeta2, k_hat, r1, r2, dl, alphaIn, Is1stElement, F_B0, F_B1, F_B2, ErrStat, ErrMsg ) + SUBROUTINE getElementHstLds_Mod1(p, m, mem, Time, pos1, pos2, Zeta1, Zeta2, k_hat, r1, r2, dl, alphaIn, Is1stElement, F_B0, F_B1, F_B2, ErrStat, ErrMsg ) + TYPE(Morison_ParameterType), INTENT( IN ) :: p + TYPE(Morison_MiscVarType), INTENT( INOUT ) :: m TYPE(Morison_MemberType), intent(in) :: mem REAL(DbKi), INTENT( IN ) :: Time REAL(ReKi), INTENT( IN ) :: pos1(3) @@ -5667,6 +5703,7 @@ SUBROUTINE getElementHstLds_Mod1( mem, Time, pos1, pos2, Zeta1, Zeta2, k_hat, r1 REAL(ReKi) :: h0, rh, a0, b0, a0b0, s0, C_1, C_2, Z0 REAL(ReKi) :: sinGamma, cosGamma, tanGamma REAL(ReKi) :: FbVec(3), MbVec(3), FSInt(3), n_hat(3), t_hat(3), s_hat(3), r_hat(3) + REAL(ReKi) :: z1, z2 INTEGER(IntKi) :: errStat2 CHARACTER(ErrMsgLen) :: errMsg2 INTEGER(IntKi), PARAMETER :: pwr = 3 ! Exponent for buoyancy node distribution smoothing @@ -5678,6 +5715,9 @@ SUBROUTINE getElementHstLds_Mod1( mem, Time, pos1, pos2, Zeta1, Zeta2, k_hat, r1 dRdl = (r2 - r1)/dl + z1 = pos1(3) + z2 = pos2(3) + IF ( (z1 < Zeta1) .AND. (z2 < Zeta2) ) THEN ! If element is fully submerged ! Compute the waterplane shape, the submerged volume, and it's geometric center ! No need to consider tapered and non-tapered elements separately @@ -5686,11 +5726,11 @@ SUBROUTINE getElementHstLds_Mod1( mem, Time, pos1, pos2, Zeta1, Zeta2, k_hat, r1 ! Hydrostatic force on element FbVec = (/0.0_ReKi,0.0_ReKi,Vs/) - Pi*( r2*r2*z2 - r1*r1*z1) *k_hat - FbVec = p%WaveField%WtrDens * g * FbVec + FbVec = p%WaveField%WtrDens * p%Gravity * FbVec ! Hydrostatic moment on element about the lower node MbVec = (Vhc+0.25*Pi*(r2**4-r1**4)) * Cross_Product(k_hat,(/0.0_ReKi,0.0_ReKi,1.0_ReKi/)) - MbVec = p%WaveField%WtrDens * g * MbVec + MbVec = p%WaveField%WtrDens * p%Gravity * MbVec ! Distribute element load to nodes alpha = alphaIn*(z2-Zeta2)**pwr/(alphaIn*(z2-Zeta2)**pwr+(1.0_ReKi-alphaIn)*(z1-Zeta1)**pwr) @@ -5715,7 +5755,7 @@ SUBROUTINE getElementHstLds_Mod1( mem, Time, pos1, pos2, Zeta1, Zeta2, k_hat, r1 rh = r1 + h0*dRdl ! Estimate the free-surface normal at the free-surface intersection, n_hat IF ( p%WaveField%WaveStMod > 0_IntKi ) THEN ! If wave stretching is enabled, compute free surface normal - CALL GetFreeSurfaceNormal( Time, FSInt, rh, n_hat, ErrStat2, ErrMsg2 ) + CALL GetFreeSurfaceNormal(p, m, Time, FSInt, rh, n_hat, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) ELSE ! Without wave stretching, use the normal of the SWL n_hat = (/0.0_ReKi,0.0_ReKi,1.0_ReKi/) @@ -5776,13 +5816,13 @@ SUBROUTINE getElementHstLds_Mod1( mem, Time, pos1, pos2, Zeta1, Zeta2, k_hat, r1 ! Hydrostatic force on element FbVec = (/0.0_ReKi,0.0_ReKi,Vs/) - Pi*a0b0*Z0*n_hat + Pi*r1**2*z1*k_hat - FbVec = p%WaveField%WtrDens * g * FbVec + FbVec = p%WaveField%WtrDens * p%Gravity * FbVec ! Hydrostatic moment on element about the lower node MbVec = Cross_Product( Vrc*r_hat+Vhc*k_hat, (/0.0_ReKi,0.0_ReKi,1.0_ReKi/) ) & + 0.25*Pi*a0b0* ( ( s_hat(3)*a0*a0 + 4.0*(s0-h0*sinGamma)*Z0 )*t_hat - t_hat(3)*b0*b0*s_hat ) & - 0.25*Pi*r1**4*( r_hat(3) *t_hat - t_hat(3) * r_hat ) - MbVec = p%WaveField%WtrDens * g * MbVec + MbVec = p%WaveField%WtrDens * p%Gravity * MbVec IF ( Is1stElement ) THEN ! This is the 1st element of the member ! Assign the element load to the lower (1st) node of the member @@ -5802,7 +5842,8 @@ SUBROUTINE getElementHstLds_Mod1( mem, Time, pos1, pos2, Zeta1, Zeta2, k_hat, r1 END IF END SUBROUTINE getElementHstLds_Mod1 - SUBROUTINE YawJoint(JointNo,PtfmRefY,AM_End,An_End,DP_Const_End,I_MG_End,ErrStat,ErrMsg) + SUBROUTINE YawJoint(p, JointNo, PtfmRefY, AM_End, An_End, DP_Const_End, I_MG_End, ErrStat, ErrMsg) + TYPE(Morison_ParameterType), INTENT( IN ) :: p Integer(IntKi), intent(in ) :: JointNo Real(ReKi), intent(in ) :: PtfmRefY Real(ReKi), intent( out) :: AM_End(3,3) @@ -5831,8 +5872,11 @@ SUBROUTINE YawJoint(JointNo,PtfmRefY,AM_End,An_End,DP_Const_End,I_MG_End,ErrStat END SUBROUTINE YawJoint - SUBROUTINE getMemBallastHiPt(member,z_hi, ErrStat, ErrMsg) + SUBROUTINE getMemBallastHiPt(p, m, u, member, z_hi, ErrStat, ErrMsg) ! This subroutine returns the highest point of a member's internal ballast + Type(Morison_ParameterType), intent(in ) :: p + Type(Morison_MiscVarType), intent(in ) :: m + Type(Morison_InputType), intent(in ) :: u Type(Morison_MemberType), intent(in ) :: member Real(ReKi), intent( out) :: z_hi Integer(IntKi), intent( out) :: ErrStat @@ -5938,8 +5982,12 @@ SUBROUTINE getMemBallastHiPt(member,z_hi, ErrStat, ErrMsg) END SUBROUTINE getMemBallastHiPt - SUBROUTINE GetDistDrag_Rec(Time,mem,i,dSadl_p,dSbdl_p,f_hydro,ErrStat,ErrMsg,SubRatio,vrelFSInt) + SUBROUTINE GetDistDrag_Rec(p, m, u, xd, Time, mem, i, dSadl_p, dSbdl_p, f_hydro, ErrStat, ErrMsg, SubRatio, vrelFSInt) ! Compute the distributed (axial and transverse) drag per unit length for rectangular sections + TYPE(Morison_ParameterType), intent(in ) :: p !< Morison parameters + Type(Morison_MiscVarType),intent(inout) :: m !< Miscellaneous variables + Type(Morison_InputType) , intent(in ) :: u !< Morison inputs + Type(Morison_DiscreteStateType), intent(in ) :: xd !< Current discrete state Real(DbKi) , intent(in ) :: Time !< Current simulation time in seconds Type(Morison_MemberType), intent(in ) :: mem !< Current member Integer(IntKi) , intent(in ) :: i !< Node number within the member (not the global node index) @@ -6107,16 +6155,6 @@ SUBROUTINE GetDistDrag_Rec(Time,mem,i,dSadl_p,dSbdl_p,f_hydro,ErrStat,ErrMsg,Sub END SUBROUTINE GetDistDrag_Rec - - logical function Failed() - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - Failed = ErrStat >= AbortErrLev - !if (Failed) then - ! call FailCleanup() - !endif - end function Failed - -END SUBROUTINE Morison_CalcOutput !---------------------------------------------------------------------------------------------------------------------------------- subroutine LumpDistrHydroLoads( f_hydro, k_hat, dl, h_c, lumpedLoad ) real(ReKi), intent(in ) :: f_hydro(3) From dade1a0e5f828d6563092ba0d985f56d8f400673 Mon Sep 17 00:00:00 2001 From: Derek Slaughter Date: Tue, 30 Dec 2025 15:55:25 -0500 Subject: [PATCH 03/12] Modified SetErrStat to overwrite contents of ErrMess if ErrStat is ErrID_None when entering the function. Therefore, each function doesn't need to set ErrMsg = "" at the beginning which reduces some overhead. --- modules/nwtc-library/src/NWTC_Base.f90 | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/modules/nwtc-library/src/NWTC_Base.f90 b/modules/nwtc-library/src/NWTC_Base.f90 index d6b45cfe48..33235ec261 100644 --- a/modules/nwtc-library/src/NWTC_Base.f90 +++ b/modules/nwtc-library/src/NWTC_Base.f90 @@ -100,9 +100,12 @@ pure subroutine SetErrStat (ErrStatLcl, ErrMessLcl, ErrStat, ErrMess, RoutineNam CHARACTER(*), INTENT(IN ) :: RoutineName ! Name of the routine error occurred in IF ( ErrStatLcl /= ErrID_None ) THEN - IF (ErrStat /= ErrID_None) ErrMess = TRIM(ErrMess)//new_line('a') - ErrMess = TRIM(ErrMess)//TRIM(RoutineName)//':'//TRIM(ErrMessLcl) - ErrStat = MAX(ErrStat,ErrStatLcl) + IF (ErrStat /= ErrID_None) then + ErrMess = TRIM(ErrMess)//new_line('a')//TRIM(RoutineName)//':'//TRIM(ErrMessLcl) + else + ErrMess = TRIM(RoutineName)//':'//TRIM(ErrMessLcl) + END IF + ErrStat = MAX(ErrStat, ErrStatLcl) END IF end subroutine From 228d47e0c56a9d8864cc397a2526f6d1ccbf4d62 Mon Sep 17 00:00:00 2001 From: Derek Slaughter Date: Tue, 30 Dec 2025 15:55:50 -0500 Subject: [PATCH 04/12] Remove some instances of ErrMsg = "" --- modules/inflowwind/src/IfW_FlowField.f90 | 2 -- .../src/NetLib/lapack/NWTC_LAPACK.f90 | 1 - modules/seastate/src/SeaSt_WaveField.f90 | 19 +++---------------- 3 files changed, 3 insertions(+), 19 deletions(-) diff --git a/modules/inflowwind/src/IfW_FlowField.f90 b/modules/inflowwind/src/IfW_FlowField.f90 index 4b04737dc1..586452aa40 100644 --- a/modules/inflowwind/src/IfW_FlowField.f90 +++ b/modules/inflowwind/src/IfW_FlowField.f90 @@ -72,7 +72,6 @@ subroutine IfW_FlowField_GetVelAcc(FF, IStart, Time, PositionXYZ, VelocityUVW, A logical :: GridExceedAllow ! is this point allowed to exceed bounds of wind grid ErrStat = ErrID_None - ErrMsg = "" ! Get number of points to evaluate NumPoints = size(PositionXYZ, dim=2) @@ -737,7 +736,6 @@ subroutine Grid3DField_GetCell(G3D, Time, Position, CalcAccel, AllowExtrap, & logical :: InGrid ErrStat = ErrID_None - ErrMsg = "" ! Initialize to no extrapolation (modified in bounds routines) AllExtrap = ExtrapNone diff --git a/modules/nwtc-library/src/NetLib/lapack/NWTC_LAPACK.f90 b/modules/nwtc-library/src/NetLib/lapack/NWTC_LAPACK.f90 index bf48d36268..2ee568ebb1 100644 --- a/modules/nwtc-library/src/NetLib/lapack/NWTC_LAPACK.f90 +++ b/modules/nwtc-library/src/NetLib/lapack/NWTC_LAPACK.f90 @@ -574,7 +574,6 @@ SUBROUTINE LAPACK_SGEMM( TRANSA, TRANSB, ALPHA, A, B, BETA, C, ErrStat, ErrMsg ) END IF ErrStat = ErrID_None - ErrMsg = "" IF ( K /= Kb ) THEN diff --git a/modules/seastate/src/SeaSt_WaveField.f90 b/modules/seastate/src/SeaSt_WaveField.f90 index 6e48c1dea3..ad42da2264 100644 --- a/modules/seastate/src/SeaSt_WaveField.f90 +++ b/modules/seastate/src/SeaSt_WaveField.f90 @@ -39,7 +39,6 @@ FUNCTION WaveField_GetNodeTotalWaveElev( WaveField, WaveField_m, Time, pos, ErrS character(ErrMsgLen) :: errMsg2 ErrStat = ErrID_None - ErrMsg = "" IF (ALLOCATED(WaveField%WaveElev1) .or. ALLOCATED(WaveField%WaveElev2)) then CALL WaveField_Interp_Setup3D(Time, pos, WaveField%GridParams, WaveField_m, ErrStat2, ErrMsg2) @@ -84,7 +83,6 @@ SUBROUTINE WaveField_GetNodeWaveNormal( WaveField, WaveField_m, Time, pos, r, n, character(ErrMsgLen) :: errMsg2 ErrStat = ErrID_None - ErrMsg = "" r1 = MAX(r,real(1.0e-6,ReKi)) ! In case r is zero @@ -134,7 +132,6 @@ SUBROUTINE WaveField_GetNodeWaveKin( WaveField, WaveField_m, Time, pos, forceNod character(ErrMsgLen) :: errMsg2 ErrStat = ErrID_None - ErrMsg = "" posXY = pos(1:2) posXY0 = (/pos(1),pos(2),0.0_ReKi/) @@ -282,7 +279,6 @@ SUBROUTINE WaveField_GetNodeWaveVel( WaveField, WaveField_m, Time, pos, forceNod character(ErrMsgLen) :: errMsg2 ErrStat = ErrID_None - ErrMsg = "" posXY = pos(1:2) posXY0 = (/pos(1),pos(2),0.0_ReKi/) @@ -396,7 +392,6 @@ SUBROUTINE WaveField_GetNodeWaveVelAcc( WaveField, WaveField_m, Time, pos, force character(ErrMsgLen) :: errMsg2 ErrStat = ErrID_None - ErrMsg = "" posXY = pos(1:2) posXY0 = (/pos(1),pos(2),0.0_ReKi/) @@ -526,7 +521,6 @@ SUBROUTINE WaveField_GetWaveKin( WaveField, WaveField_m, Time, pos, forceNodeInW real(ReKi), allocatable :: FV_DC(:,:), FA_DC(:,:) ErrStat = ErrID_None - ErrMsg = "" NumPoints = size(pos, dim=2) DO i = 1, NumPoints @@ -561,9 +555,9 @@ logical function Failed() call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) Failed = ErrStat >= AbortErrLev end function - logical function FailedMsg(ErrMsg2) - character(*), intent(in ) :: ErrMsg2 - call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + logical function FailedMsg(ErrMsgTmp) + character(*), intent(in ) :: ErrMsgTmp + call SetErrStat( ErrStat2, ErrMsgTmp, ErrStat, ErrMsg, RoutineName ) FailedMsg = ErrStat >= AbortErrLev end function end subroutine WaveField_GetWaveKin @@ -591,7 +585,6 @@ SUBROUTINE WaveField_GetWaveVelAcc_AD( WaveField, WaveField_m, StartNode, Time, character(ErrMsgLen) :: errMsg2 ErrStat = ErrID_None - ErrMsg = "" MSL2SWL = WaveField%MSL2SWL WtrDpth = WaveField%EffWtrDpth - MSL2SWL @@ -668,7 +661,6 @@ SUBROUTINE WaveField_GetMeanDynSurfCurr( WaveField, WaveTMax, WaveDT, CurrVxi0, character(ErrMsgLen) :: errMsg2 ErrStat = ErrID_None - ErrMsg = "" CurrVxi0 = 0.0_SiKi CurrVyi0 = 0.0_SiKi @@ -726,7 +718,6 @@ subroutine SetCartesianXYIndex(p, pZero, delta, nMax, Indx_Lo, Indx_Hi, isopc, F real(ReKi) :: Tmp ErrStat = ErrID_None - ErrMsg = "" isopc = -1.0 Indx_Lo = 0 @@ -792,7 +783,6 @@ subroutine SetCartesianZIndex(p, z_depth, delta, nMax, Indx_Lo, Indx_Hi, isopc, real(ReKi) :: Tmp ErrStat = ErrID_None - ErrMsg = "" isopc = -1.0 Indx_Lo = 0 @@ -849,7 +839,6 @@ subroutine SetTimeIndex(Time, deltaT, nMax, Indx_Lo, Indx_Hi, isopc, ErrStat, Er real(ReKi) :: Tmp ErrStat = ErrID_None - ErrMsg = "" isopc = -1.0 Indx_Lo = 0 @@ -905,7 +894,6 @@ subroutine WaveField_Interp_Setup4D( Time, Position, p, m, ErrStat, ErrMsg ) character(ErrMsgLen) :: ErrMsg2 ErrStat = ErrID_None - ErrMsg = "" ! Find the bounding indices for time call SetTimeIndex(Time, p%delta(1), p%n(1), m%Indx_Lo(1), m%Indx_Hi(1), isopc(1), ErrStat2, ErrMsg2) @@ -974,7 +962,6 @@ subroutine WaveField_Interp_Setup3D( Time, Position, p, m, ErrStat, ErrMsg ) character(ErrMsgLen) :: ErrMsg2 ErrStat = ErrID_None - ErrMsg = "" ! Find the bounding indices for time call SetTimeIndex(Time, p%delta(1), p%n(1), m%Indx_Lo(1), m%Indx_Hi(1), isopc(1), ErrStat2, ErrMsg2) From 372149998ea960bf71e100984d7fb4a4f736059f Mon Sep 17 00:00:00 2001 From: Derek Slaughter Date: Tue, 30 Dec 2025 15:56:13 -0500 Subject: [PATCH 05/12] Restructure getWaterKin in attempt to be more performant --- modules/moordyn/src/MoorDyn_Misc.f90 | 111 +++++++++++++-------------- 1 file changed, 52 insertions(+), 59 deletions(-) diff --git a/modules/moordyn/src/MoorDyn_Misc.f90 b/modules/moordyn/src/MoorDyn_Misc.f90 index 65fd5edddd..56ae66cf2d 100644 --- a/modules/moordyn/src/MoorDyn_Misc.f90 +++ b/modules/moordyn/src/MoorDyn_Misc.f90 @@ -948,62 +948,42 @@ SUBROUTINE getWaterKin(p, m, x, y, z, t, U, Ud, zeta, PDyn, ErrStat, ErrMsg) ErrStat = ErrID_None ErrMsg = "" - IF (p%WaterKin == 3 .AND. (.NOT. m%IC_gen)) THEN ! disable wavekin 3 during IC_gen, otherwise will never find steady state (becasue of waves) - - ! SeaState throws warning when queried location is out of bounds from the SeaState grid, so no need to handle here + ! Initialize values to zero + U = 0.0_DbKi + Ud = 0.0_DbKi + zeta = 0.0_DbKi + PDyn = 0.0_DbKi - ! Pack all MD inputs to WaveGrid input data types (double to single) - ! (only pos needed becasue time is double in wave field, all other are outputs that will be set by WaveField_GetNodeWaveKin) - xyz_sp = REAL((/ x, y, z /),SiKi) + select case (p%WaterKin) - ! for now we will force the node to be in the water (forceNodeInWater = True). Rods handle partial submergence seperately so they need to get information from SeaState - CALL WaveField_GetNodeWaveKin(p%WaveField, m%WaveField_m, t, xyz_sp, .TRUE., .TRUE., nodeInWater, WaveElev1, WaveElev2, zeta_sp, PDyn_sp, U_sp, Ud_sp, FAMCF, ErrStat2, ErrMsg2 ) ! outputs: nodeInWater, WaveElev1, WaveElev2, FAMCF all unused - CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + ! old or hybrid approach. SeaState contributions handeled in setupWaterKin, just proceed using old method + case (1,2) - ! Unpack all WaveGrid outputs to MD output data types (single to double) - U = REAL(U_sp,DbKi) - Ud = REAL(Ud_sp,DbKi) - zeta = REAL(zeta_sp,DbKi) - PDyn = REAL(PDyn_sp,DbKi) - - ELSEIF (p%WaterKin == 1 .OR. p%WaterKin == 2) THEN ! old or hybrid approach. SeaState contributions handeled in setupWaterKin, just proceed using old method - - ! If wave kinematics enabled, get interpolated values from grid - IF (p%WaveKin > 0) THEN - - ! find time interpolation indices and coefficients - it = floor(t/ p%dtWave) + 1 ! add 1 because Fortran indexing starts at 1 - ft = (t - (it-1)*p%dtWave)/p%dtWave - m%WaveTi = it + ! find time interpolation indices and coefficients + it = floor(t/ p%dtWave) + 1 ! add 1 because Fortran indexing starts at 1 + ft = (t - (it-1)*p%dtWave)/p%dtWave + m%WaveTi = it + + ! find x-y interpolation indices and coefficients + CALL getInterpNumsSiKi(p%pxWave , REAL(x,SiKi), 1, ix, fx) ! wave grid + CALL getInterpNumsSiKi(p%pyWave , REAL(y,SiKi), 1, iy, fy) ! wave grid - ! find x-y interpolation indices and coefficients - CALL getInterpNumsSiKi(p%pxWave , REAL(x,SiKi), 1, ix, fx) ! wave grid - CALL getInterpNumsSiKi(p%pyWave , REAL(y,SiKi), 1, iy, fy) ! wave grid - - ! interpolate wave elevation - CALL calculate3Dinterpolation(p%zeta, ix, iy, it, fx, fy, ft, zeta) - - ! compute modified z coordinate to be used for interpolating velocities and accelerations with Wheeler stretching - zp = ( z - zeta ) * p%WtrDpth/( p%WtrDpth + zeta ) - - CALL getInterpNumsSiKi(p%pzWave , REAL(zp,SiKi), 1, iz, fz) ! wave grid - - ! interpolate everything else - CALL calculate4Dinterpolation(p%PDyn , ix, iy, iz, it, fx, fy, fz, ft, PDyn) - CALL calculate4Dinterpolation(p%uxWave, ix, iy, iz, it, fx, fy, fz, ft, U(1) ) - CALL calculate4Dinterpolation(p%uyWave, ix, iy, iz, it, fx, fy, fz, ft, U(2) ) - CALL calculate4Dinterpolation(p%uzWave, ix, iy, iz, it, fx, fy, fz, ft, U(3) ) - CALL calculate4Dinterpolation(p%axWave, ix, iy, iz, it, fx, fy, fz, ft, Ud(1) ) - CALL calculate4Dinterpolation(p%ayWave, ix, iy, iz, it, fx, fy, fz, ft, Ud(2) ) - CALL calculate4Dinterpolation(p%azWave, ix, iy, iz, it, fx, fy, fz, ft, Ud(3) ) + ! interpolate wave elevation + CALL calculate3Dinterpolation(p%zeta, ix, iy, it, fx, fy, ft, zeta) + + ! compute modified z coordinate to be used for interpolating velocities and accelerations with Wheeler stretching + zp = ( z - zeta ) * p%WtrDpth/( p%WtrDpth + zeta ) - ELSE ! set things to zero if wave kinematics not enabled - U = 0.0_DbKi - Ud = 0.0_DbKi - zeta = 0.0_DbKi - PDyn = 0.0_DbKi + CALL getInterpNumsSiKi(p%pzWave , REAL(zp,SiKi), 1, iz, fz) ! wave grid - ENDIF + ! interpolate everything else + CALL calculate4Dinterpolation(p%PDyn , ix, iy, iz, it, fx, fy, fz, ft, PDyn) + CALL calculate4Dinterpolation(p%uxWave, ix, iy, iz, it, fx, fy, fz, ft, U(1) ) + CALL calculate4Dinterpolation(p%uyWave, ix, iy, iz, it, fx, fy, fz, ft, U(2) ) + CALL calculate4Dinterpolation(p%uzWave, ix, iy, iz, it, fx, fy, fz, ft, U(3) ) + CALL calculate4Dinterpolation(p%axWave, ix, iy, iz, it, fx, fy, fz, ft, Ud(1) ) + CALL calculate4Dinterpolation(p%ayWave, ix, iy, iz, it, fx, fy, fz, ft, Ud(2) ) + CALL calculate4Dinterpolation(p%azWave, ix, iy, iz, it, fx, fy, fz, ft, Ud(3) ) ! If current kinematics enabled, add interpolated current values from profile IF (p%Current > 0) THEN @@ -1021,15 +1001,28 @@ SUBROUTINE getWaterKin(p, m, x, y, z, t, U, Ud, zeta, PDyn, ErrStat, ErrMsg) U(2) = U(2) + (1.0-fz)*p%uyCurrent(iz0) + fz*p%uyCurrent(iz1) END IF - ELSEIF (p%WaterKin > 3) THEN - CALL SetErrStat(ErrID_Fatal, "Invalid p%WaterKin value found in getWaterKin", ErrStat, ErrMsg, RoutineName) - - ELSE ! set things to zero if Water Kinematics not enabled - U = 0.0_DbKi - Ud = 0.0_DbKi - zeta = 0.0_DbKi - PDyn = 0.0_DbKi - ENDIF + case (3) + + ! Disable wavekin 3 during IC_gen, otherwise will never find steady state (becasue of waves) + if (m%IC_gen) return + + ! SeaState throws warning when queried location is out of bounds from the SeaState grid, so no need to handle here + + ! Pack all MD inputs to WaveGrid input data types (double to single) + ! (only pos needed becasue time is double in wave field, all other are outputs that will be set by WaveField_GetNodeWaveKin) + xyz_sp = REAL((/ x, y, z /),SiKi) + + ! for now we will force the node to be in the water (forceNodeInWater = True). Rods handle partial submergence seperately so they need to get information from SeaState + CALL WaveField_GetNodeWaveKin(p%WaveField, m%WaveField_m, t, xyz_sp, .TRUE., .TRUE., nodeInWater, WaveElev1, WaveElev2, zeta_sp, PDyn_sp, U_sp, Ud_sp, FAMCF, ErrStat2, ErrMsg2 ) ! outputs: nodeInWater, WaveElev1, WaveElev2, FAMCF all unused + CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + + ! Unpack all WaveGrid outputs to MD output data types (single to double) + U = REAL(U_sp,DbKi) + Ud = REAL(Ud_sp,DbKi) + zeta = REAL(zeta_sp,DbKi) + PDyn = REAL(PDyn_sp,DbKi) + + end select END SUBROUTINE getWaterKin From 682224bcefe82ad36e25e760ce8c84b27dd7e79c Mon Sep 17 00:00:00 2001 From: Derek Slaughter Date: Wed, 31 Dec 2025 16:30:51 +0000 Subject: [PATCH 06/12] Fix bad refactored logic in getWaterKin --- modules/moordyn/src/MoorDyn_Misc.f90 | 81 +++++++++++++++------------- 1 file changed, 45 insertions(+), 36 deletions(-) diff --git a/modules/moordyn/src/MoorDyn_Misc.f90 b/modules/moordyn/src/MoorDyn_Misc.f90 index 56ae66cf2d..03aed01572 100644 --- a/modules/moordyn/src/MoorDyn_Misc.f90 +++ b/modules/moordyn/src/MoorDyn_Misc.f90 @@ -948,72 +948,79 @@ SUBROUTINE getWaterKin(p, m, x, y, z, t, U, Ud, zeta, PDyn, ErrStat, ErrMsg) ErrStat = ErrID_None ErrMsg = "" - ! Initialize values to zero + ! Initialize outputs to zero U = 0.0_DbKi Ud = 0.0_DbKi zeta = 0.0_DbKi PDyn = 0.0_DbKi select case (p%WaterKin) + ! no wave kinematics, do nothing (already initialized to zero above) + case (0) ! old or hybrid approach. SeaState contributions handeled in setupWaterKin, just proceed using old method case (1,2) - ! find time interpolation indices and coefficients - it = floor(t/ p%dtWave) + 1 ! add 1 because Fortran indexing starts at 1 - ft = (t - (it-1)*p%dtWave)/p%dtWave - m%WaveTi = it - - ! find x-y interpolation indices and coefficients - CALL getInterpNumsSiKi(p%pxWave , REAL(x,SiKi), 1, ix, fx) ! wave grid - CALL getInterpNumsSiKi(p%pyWave , REAL(y,SiKi), 1, iy, fy) ! wave grid - - ! interpolate wave elevation - CALL calculate3Dinterpolation(p%zeta, ix, iy, it, fx, fy, ft, zeta) - - ! compute modified z coordinate to be used for interpolating velocities and accelerations with Wheeler stretching - zp = ( z - zeta ) * p%WtrDpth/( p%WtrDpth + zeta ) - - CALL getInterpNumsSiKi(p%pzWave , REAL(zp,SiKi), 1, iz, fz) ! wave grid + ! If wave kinematics enabled, get interpolated values from grid + if (p%WaveKin > 0) then - ! interpolate everything else - CALL calculate4Dinterpolation(p%PDyn , ix, iy, iz, it, fx, fy, fz, ft, PDyn) - CALL calculate4Dinterpolation(p%uxWave, ix, iy, iz, it, fx, fy, fz, ft, U(1) ) - CALL calculate4Dinterpolation(p%uyWave, ix, iy, iz, it, fx, fy, fz, ft, U(2) ) - CALL calculate4Dinterpolation(p%uzWave, ix, iy, iz, it, fx, fy, fz, ft, U(3) ) - CALL calculate4Dinterpolation(p%axWave, ix, iy, iz, it, fx, fy, fz, ft, Ud(1) ) - CALL calculate4Dinterpolation(p%ayWave, ix, iy, iz, it, fx, fy, fz, ft, Ud(2) ) - CALL calculate4Dinterpolation(p%azWave, ix, iy, iz, it, fx, fy, fz, ft, Ud(3) ) + ! find time interpolation indices and coefficients + it = floor(t/ p%dtWave) + 1 ! add 1 because Fortran indexing starts at 1 + ft = (t - (it-1)*p%dtWave)/p%dtWave + m%WaveTi = it + + ! find x-y interpolation indices and coefficients + CALL getInterpNumsSiKi(p%pxWave, REAL(x, SiKi), 1, ix, fx) ! wave grid + CALL getInterpNumsSiKi(p%pyWave, REAL(y, SiKi), 1, iy, fy) ! wave grid + + ! interpolate wave elevation + CALL calculate3Dinterpolation(p%zeta, ix, iy, it, fx, fy, ft, zeta) + + ! compute modified z coordinate to be used for interpolating velocities and accelerations with Wheeler stretching + zp = (z - zeta) * p%WtrDpth/(p%WtrDpth + zeta) + + CALL getInterpNumsSiKi(p%pzWave, REAL(zp, SiKi), 1, iz, fz) ! wave grid + + ! interpolate everything else + CALL calculate4Dinterpolation(p%PDyn , ix, iy, iz, it, fx, fy, fz, ft, PDyn) + CALL calculate4Dinterpolation(p%uxWave, ix, iy, iz, it, fx, fy, fz, ft, U(1)) + CALL calculate4Dinterpolation(p%uyWave, ix, iy, iz, it, fx, fy, fz, ft, U(2)) + CALL calculate4Dinterpolation(p%uzWave, ix, iy, iz, it, fx, fy, fz, ft, U(3)) + CALL calculate4Dinterpolation(p%axWave, ix, iy, iz, it, fx, fy, fz, ft, Ud(1)) + CALL calculate4Dinterpolation(p%ayWave, ix, iy, iz, it, fx, fy, fz, ft, Ud(2)) + CALL calculate4Dinterpolation(p%azWave, ix, iy, iz, it, fx, fy, fz, ft, Ud(3)) + end if ! If current kinematics enabled, add interpolated current values from profile - IF (p%Current > 0) THEN + if (p%Current > 0) then - CALL getInterpNumsSiKi(p%pzCurrent, REAL(z,SiKi), 1, iz0, fz) + CALL getInterpNumsSiKi(p%pzCurrent, REAL(z, SiKi), 1, iz0, fz) - IF (fz == 0) THEN ! handle end case conditions + if (fz == 0) then ! handle end case conditions iz1 = iz0 - ELSE - iz1 = min(iz0+1,size(p%pzCurrent)) ! don't overstep bounds - END IF + else + iz1 = min(iz0+1, size(p%pzCurrent)) ! don't overstep bounds + end if ! Add the current velocities to the wave velocities (if any) U(1) = U(1) + (1.0-fz)*p%uxCurrent(iz0) + fz*p%uxCurrent(iz1) U(2) = U(2) + (1.0-fz)*p%uyCurrent(iz0) + fz*p%uyCurrent(iz1) - END IF + end if + ! SeaState wave kinematics case (3) - ! Disable wavekin 3 during IC_gen, otherwise will never find steady state (becasue of waves) + ! disable wavekin 3 during IC_gen, otherwise will never find steady state (becasue of waves) if (m%IC_gen) return ! SeaState throws warning when queried location is out of bounds from the SeaState grid, so no need to handle here ! Pack all MD inputs to WaveGrid input data types (double to single) - ! (only pos needed becasue time is double in wave field, all other are outputs that will be set by WaveField_GetNodeWaveKin) - xyz_sp = REAL((/ x, y, z /),SiKi) + ! (only pos needed because time is double in wave field, all other are outputs that will be set by WaveField_GetNodeWaveKin) + xyz_sp = REAL([x, y, z], SiKi) ! for now we will force the node to be in the water (forceNodeInWater = True). Rods handle partial submergence seperately so they need to get information from SeaState - CALL WaveField_GetNodeWaveKin(p%WaveField, m%WaveField_m, t, xyz_sp, .TRUE., .TRUE., nodeInWater, WaveElev1, WaveElev2, zeta_sp, PDyn_sp, U_sp, Ud_sp, FAMCF, ErrStat2, ErrMsg2 ) ! outputs: nodeInWater, WaveElev1, WaveElev2, FAMCF all unused + CALL WaveField_GetNodeWaveKin(p%WaveField, m%WaveField_m, t, xyz_sp, .true., .true., nodeInWater, WaveElev1, WaveElev2, zeta_sp, PDyn_sp, U_sp, Ud_sp, FAMCF, ErrStat2, ErrMsg2 ) ! outputs: nodeInWater, WaveElev1, WaveElev2, FAMCF all unused CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) ! Unpack all WaveGrid outputs to MD output data types (single to double) @@ -1022,6 +1029,8 @@ SUBROUTINE getWaterKin(p, m, x, y, z, t, U, Ud, zeta, PDyn, ErrStat, ErrMsg) zeta = REAL(zeta_sp,DbKi) PDyn = REAL(PDyn_sp,DbKi) + case default + call SetErrStat(ErrID_Fatal, "Invalid value of p%WaterKin", ErrStat, ErrMsg, RoutineName) end select END SUBROUTINE getWaterKin From dea020b20b1a88e07317bf6945a84628883d11af Mon Sep 17 00:00:00 2001 From: Derek Slaughter Date: Wed, 31 Dec 2025 12:18:42 -0500 Subject: [PATCH 07/12] Minor optimization in FVW ui_seg (~6%) --- modules/aerodyn/src/FVW_BiotSavart.f90 | 96 ++++++++++++++------------ 1 file changed, 51 insertions(+), 45 deletions(-) diff --git a/modules/aerodyn/src/FVW_BiotSavart.f90 b/modules/aerodyn/src/FVW_BiotSavart.f90 index 36c73d1dd9..038a38a953 100644 --- a/modules/aerodyn/src/FVW_BiotSavart.f90 +++ b/modules/aerodyn/src/FVW_BiotSavart.f90 @@ -134,6 +134,7 @@ subroutine ui_seg(iCPStart, iCPEnd, CPs, & real(ReKi) :: norm2_r0 !< Squared length of the segment d = (r1xr2)/r0 real(ReKi) :: xa, ya, za, xb, yb, zb !< Coordinates of X-Xa and X-Xb real(ReKi) :: exp_value !< + real(ReKi) :: CPs_icp(3) !< ! Branching based on regularization model ! NOTE: copy paste of code is done for optimization! @@ -141,14 +142,15 @@ subroutine ui_seg(iCPStart, iCPEnd, CPs, & select case (RegFunction) case ( idRegNone ) ! No vortex core !$OMP PARALLEL default(shared) - !$OMP do private(icp,is,Uind,P1,P2,crossprod,denominator,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb) schedule(runtime) + !$OMP do private(icp,is,CPs_icp,Uind,P1,P2,crossprod,denominator,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb) schedule(runtime) do icp=iCPStart,iCPEnd ! loop on CPs + Uind = 0.0_ReKi + CPs_icp = CPs(:,icp) do is=iSegStart,iSegEnd ! loop on selected segments - Uind = 0.0_ReKi P1 = SegPoints(1:3, SegConnct(1,is)) ! Segment extremity points P2 = SegPoints(1:3, SegConnct(2,is)) - xa=CPs(1,icp)-P1(1); ya=CPs(2,icp)-P1(2); za=CPs(3,icp)-P1(3); - xb=CPs(1,icp)-P2(1); yb=CPs(2,icp)-P2(2); zb=CPs(3,icp)-P2(3); + xa=CPs_icp(1)-P1(1); ya=CPs_icp(2)-P1(2); za=CPs_icp(3)-P1(3); + xb=CPs_icp(1)-P2(1); yb=CPs_icp(2)-P2(2); zb=CPs_icp(3)-P2(3); norm_a = sqrt(xa*xa + ya*ya + za*za) norm_b = sqrt(xb*xb + yb*yb + zb*zb) denominator = norm_a*norm_b*(norm_a*norm_b + xa*xb+ya*yb+za*zb) @@ -162,26 +164,27 @@ subroutine ui_seg(iCPStart, iCPEnd, CPs, & ! --- Far field TODO ! --- NO Regularization (close field) Kv = SegGamma(is)*fourpi_inv*(norm_a+norm_b)/(denominator + MINDENOM) - Uind(1:3) = Kv*crossprod(1:3) + Uind(1:3) = Uind(1:3) + Kv*crossprod(1:3) end if end if end if - Uind_out(1:3,icp) = Uind_out(1:3,icp)+Uind(1:3) end do ! Loop on segments + Uind_out(1:3,icp) = Uind_out(1:3,icp)+Uind(1:3) enddo ! Loop on control points !$OMP END DO !$OMP END PARALLEL case ( idRegRankine ) ! Rankine !$OMP PARALLEL default(shared) - !$OMP do private(icp,is,Uind,P1,P2,crossprod,denominator,r_bar2,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb) schedule(runtime) + !$OMP do private(icp,is,CPs_icp,Uind,P1,P2,crossprod,denominator,r_bar2,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb) schedule(runtime) do icp=iCPStart,iCPEnd ! loop on CPs + Uind = 0.0_ReKi + CPs_icp = CPs(:,icp) do is=iSegStart,iSegEnd ! loop on selected segments - Uind = 0.0_ReKi P1 = SegPoints(1:3, SegConnct(1,is)) ! Segment extremity points P2 = SegPoints(1:3, SegConnct(2,is)) - xa=CPs(1,icp)-P1(1); ya=CPs(2,icp)-P1(2); za=CPs(3,icp)-P1(3); - xb=CPs(1,icp)-P2(1); yb=CPs(2,icp)-P2(2); zb=CPs(3,icp)-P2(3); + xa=CPs_icp(1)-P1(1); ya=CPs_icp(2)-P1(2); za=CPs_icp(3)-P1(3); + xb=CPs_icp(1)-P2(1); yb=CPs_icp(2)-P2(2); zb=CPs_icp(3)-P2(3); norm_a = sqrt(xa*xa + ya*ya + za*za) norm_b = sqrt(xb*xb + yb*yb + zb*zb) denominator = norm_a*norm_b*(norm_a*norm_b + xa*xb+ya*yb+za*zb) @@ -201,26 +204,27 @@ subroutine ui_seg(iCPStart, iCPEnd, CPs, & Kv=1.0_ReKi end if Kv = SegGamma(is)*fourpi_inv*Kv*(norm_a+norm_b)/(denominator + MINDENOM) - Uind(1:3) = Kv*crossprod(1:3) + Uind(1:3) = Uind(1:3) + Kv*crossprod(1:3) end if end if ! denominator size or distances too small end if ! - Uind_out(1:3,icp) = Uind_out(1:3,icp)+Uind(1:3) end do ! Loop on segments + Uind_out(1:3,icp) = Uind_out(1:3,icp) + Uind(1:3) enddo ! Loop on control points !$OMP END DO !$OMP END PARALLEL case ( idRegLambOseen ) ! LambOseen !$OMP PARALLEL default(shared) - !$OMP do private(icp,is,Uind,P1,P2,crossprod,denominator,r_bar2,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb,exp_value) schedule(runtime) + !$OMP do private(icp,is,CPs_icp,Uind,P1,P2,crossprod,denominator,r_bar2,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb,exp_value) schedule(runtime) do icp=iCPStart,iCPEnd ! loop on CPs + Uind = 0.0_ReKi + CPs_icp = CPs(:,icp) do is=iSegStart,iSegEnd ! loop on selected segments - Uind = 0.0_ReKi P1 = SegPoints(1:3, SegConnct(1,is)) ! Segment extremity points P2 = SegPoints(1:3, SegConnct(2,is)) - xa=CPs(1,icp)-P1(1); ya=CPs(2,icp)-P1(2); za=CPs(3,icp)-P1(3); - xb=CPs(1,icp)-P2(1); yb=CPs(2,icp)-P2(2); zb=CPs(3,icp)-P2(3); + xa=CPs_icp(1)-P1(1); ya=CPs_icp(2)-P1(2); za=CPs_icp(3)-P1(3); + xb=CPs_icp(1)-P2(1); yb=CPs_icp(2)-P2(2); zb=CPs_icp(3)-P2(3); norm_a = sqrt(xa*xa + ya*ya + za*za) norm_b = sqrt(xb*xb + yb*yb + zb*zb) denominator = norm_a*norm_b*(norm_a*norm_b + xa*xb+ya*yb+za*zb) @@ -241,61 +245,63 @@ subroutine ui_seg(iCPStart, iCPEnd, CPs, & Kv = 1.0_ReKi-exp(exp_value) endif Kv = SegGamma(is)*fourpi_inv*Kv*(norm_a+norm_b)/(denominator + MINDENOM) - Uind(1:3) = Kv*crossprod(1:3) + Uind(1:3) = Uind(1:3) + Kv*crossprod(1:3) endif end if ! denominator size or distances too small end if ! - Uind_out(1:3,icp) = Uind_out(1:3,icp)+Uind(1:3) end do ! Loop on segments + Uind_out(1:3,icp) = Uind_out(1:3,icp) + Uind(1:3) enddo ! Loop on control points !$OMP END DO !$OMP END PARALLEL case ( idRegVatistas ) ! Vatistas n=2 !$OMP PARALLEL default(shared) - !$OMP do private(icp,is,Uind,P1,P2,crossprod,denominator,r_bar2,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb) schedule(runtime) + !$OMP do private(icp,is,CPs_icp,Uind,P1,P2,crossprod,denominator,r_bar2,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb) schedule(runtime) do icp=iCPStart,iCPEnd ! loop on CPs + Uind = 0.0_ReKi + CPs_icp = CPs(:,icp) do is=iSegStart,iSegEnd ! loop on selected segments - Uind = 0.0_ReKi P1 = SegPoints(1:3, SegConnct(1,is)) ! Segment extremity points P2 = SegPoints(1:3, SegConnct(2,is)) - xa=CPs(1,icp)-P1(1); ya=CPs(2,icp)-P1(2); za=CPs(3,icp)-P1(3); - xb=CPs(1,icp)-P2(1); yb=CPs(2,icp)-P2(2); zb=CPs(3,icp)-P2(3); + xa=CPs_icp(1)-P1(1); ya=CPs_icp(2)-P1(2); za=CPs_icp(3)-P1(3); + xb=CPs_icp(1)-P2(1); yb=CPs_icp(2)-P2(2); zb=CPs_icp(3)-P2(3); norm_a = sqrt(xa*xa + ya*ya + za*za) norm_b = sqrt(xb*xb + yb*yb + zb*zb) denominator = norm_a*norm_b*(norm_a*norm_b + xa*xb+ya*yb+za*zb) - if (denominator>PRECISION_UI) then - crossprod(1) = ya*zb-za*yb; crossprod(2) = za*xb-xa*zb; crossprod(3) = xa*yb-ya*xb - norm2_orth = crossprod(1)**2 + crossprod(2)**2 + crossprod(3)**2 - if (norm2_orth>PRECISION_UI) then ! On the singularity, Uind(1:3)=0.0_ReKi - norm2_r0 = (xa-xb)*(xa-xb) + (ya-yb)*(ya-yb) +(za-zb)*(za-zb) - if (norm2_r0>PRECISION_UI) then - ! --- Far field TODO - ! --- Regularization (close field) --- Vatistas - norm2_orth = norm2_orth/norm2_r0 ! d = (r1xr2)/r0 - r_bar2 = norm2_orth/RegParam(is)**2 - Kv = r_bar2/sqrt(1+r_bar2**2) - Kv = SegGamma(is)*fourpi_inv*Kv*(norm_a+norm_b)/(denominator + MINDENOM) - Uind(1:3) = Kv*crossprod(1:3) - end if - end if ! denominator size or distances too small - end if ! - Uind_out(1:3,icp) = Uind_out(1:3,icp)+Uind(1:3) + ! denominator size or distances too small + if (denominator <= PRECISION_UI) cycle + crossprod(1) = ya*zb-za*yb; crossprod(2) = za*xb-xa*zb; crossprod(3) = xa*yb-ya*xb + norm2_orth = crossprod(1)**2 + crossprod(2)**2 + crossprod(3)**2 + ! On the singularity, cycle + if (norm2_orth <= PRECISION_UI) cycle + norm2_r0 = (xa-xb)*(xa-xb) + (ya-yb)*(ya-yb) +(za-zb)*(za-zb) + ! segment of zero length + if (norm2_r0 <= PRECISION_UI) cycle + ! --- Far field TODO + ! --- Regularization (close field) --- Vatistas + norm2_orth = norm2_orth/norm2_r0 ! d = (r1xr2)/r0 + r_bar2 = norm2_orth/RegParam(is)**2 + Kv = r_bar2/sqrt(1.0_ReKi+r_bar2**2) + Kv = SegGamma(is)*fourpi_inv*Kv*(norm_a+norm_b)/(denominator + MINDENOM) + Uind(1:3) = Uind(1:3) + Kv*crossprod(1:3) end do ! Loop on segments + Uind_out(1:3,icp) = Uind_out(1:3,icp) + Uind(1:3) enddo ! Loop on control points !$OMP END DO !$OMP END PARALLEL case ( idRegOffset ) ! Denominator offset !$OMP PARALLEL default(shared) - !$OMP do private(icp,is,Uind,P1,P2,crossprod,denominator,r_bar2,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb) schedule(runtime) + !$OMP do private(icp,is,CPs_icp,Uind,P1,P2,crossprod,denominator,r_bar2,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb) schedule(runtime) do icp=iCPStart,iCPEnd ! loop on CPs + Uind = 0.0_ReKi + CPs_icp = CPs(:,icp) do is=iSegStart,iSegEnd ! loop on selected segments - Uind = 0.0_ReKi P1 = SegPoints(1:3, SegConnct(1,is)) ! Segment extremity points P2 = SegPoints(1:3, SegConnct(2,is)) - xa=CPs(1,icp)-P1(1); ya=CPs(2,icp)-P1(2); za=CPs(3,icp)-P1(3); - xb=CPs(1,icp)-P2(1); yb=CPs(2,icp)-P2(2); zb=CPs(3,icp)-P2(3); + xa=CPs_icp(1)-P1(1); ya=CPs_icp(2)-P1(2); za=CPs_icp(3)-P1(3); + xb=CPs_icp(1)-P2(1); yb=CPs_icp(2)-P2(2); zb=CPs_icp(3)-P2(3); norm_a = sqrt(xa*xa + ya*ya + za*za) norm_b = sqrt(xb*xb + yb*yb + zb*zb) denominator = norm_a*norm_b*(norm_a*norm_b + xa*xb+ya*yb+za*zb) @@ -309,12 +315,12 @@ subroutine ui_seg(iCPStart, iCPEnd, CPs, & ! --- Regularization (close field) -- Offset denominator = denominator+RegParam(is)**2*norm2_r0 Kv = SegGamma(is)*fourpi_inv*(norm_a+norm_b)/(denominator + MINDENOM) - Uind(1:3) = Kv*crossprod(1:3) + Uind(1:3) = Uind(1:3) + Kv*crossprod(1:3) end if end if ! denominator size or distances too small end if ! - Uind_out(1:3,icp) = Uind_out(1:3,icp)+Uind(1:3) end do ! Loop on segments + Uind_out(1:3,icp) = Uind_out(1:3,icp)+Uind(1:3) enddo ! Loop on control points !$OMP END DO !$OMP END PARALLEL From eddeb743dc269621c805640b04f767b93bd782dd Mon Sep 17 00:00:00 2001 From: Derek Slaughter Date: Wed, 31 Dec 2025 14:49:45 -0500 Subject: [PATCH 08/12] Remove recursion from Morison RefineElementHstLds_Cyl --- modules/hydrodyn/src/Morison.f90 | 182 ++++++++++++++++++++++++------- 1 file changed, 141 insertions(+), 41 deletions(-) diff --git a/modules/hydrodyn/src/Morison.f90 b/modules/hydrodyn/src/Morison.f90 index 84b51b4ff4..a28cc215e7 100644 --- a/modules/hydrodyn/src/Morison.f90 +++ b/modules/hydrodyn/src/Morison.f90 @@ -5199,7 +5199,7 @@ SUBROUTINE getElementHstLds_Mod2_Rec(p, pos1In, pos2In, FSPtIn, k_hatIn, x_hatIn END SUBROUTINE getElementHstLds_Mod2_Rec - RECURSIVE SUBROUTINE RefineElementHstLds_Cyl(p, origin, pos1, posMid, pos2, FSPt, r1, rMid, r2, dl, dRdl,secStat1,secStatMid,secStat2, k_hat, y_hat, z_hat, n_hat, dFdl1, dFdlMid, dFdl2, recurLvl, F_B_5pt, ErrStat) + SUBROUTINE RefineElementHstLds_Cyl(p, origin, pos1, posMid, pos2, FSPt, r1, rMid, r2, dl, dRdl,secStat1,secStatMid,secStat2, k_hat, y_hat, z_hat, n_hat, dFdl1, dFdlMid, dFdl2, recurLvl, F_B_5pt, ErrStat) TYPE(Morison_ParameterType), INTENT( IN ) :: p REAL(DbKi), INTENT( IN ) :: origin(3) REAL(DbKi), INTENT( IN ) :: pos1(3) @@ -5229,63 +5229,163 @@ RECURSIVE SUBROUTINE RefineElementHstLds_Cyl(p, origin, pos1, posMid, pos2, FSPt REAL(DbKi) :: theta1,theta2 REAL(DbKi) :: posMidL(3), posMidR(3) REAL(DbKi) :: rMidL, rMidR - REAL(DbKi) :: dFdlMidL(6), dFdlMidR(6), F_B_3pt(6) + REAL(DbKi) :: dFdlMidL(6), dFdlMidR(6) + REAL(DbKi) :: F_B_3pt_sub(6), F_B_5pt_sub(6) REAL(DbKi) :: error(6), tmp(6) LOGICAL :: tolMet - INTEGER(IntKi) :: i + INTEGER(IntKi) :: i, ib INTEGER(IntKi) :: secStatMidL, secStatMidR + REAL(ReKi) :: k_hat_Re(3) + REAL(ReKi) :: y_hat_Re(3) + REAL(ReKi) :: z_hat_Re(3) + REAL(ReKi) :: n_hat_Re(3) REAL(DbKi), PARAMETER :: RelTol = 1.0E-6 REAL(DbKi), PARAMETER :: AbsTol = 1.0E-8 INTEGER(IntKi), PARAMETER :: maxRecurLvl = 50 CHARACTER(*), PARAMETER :: RoutineName = "RefineElementHstLds_Cyl" - + + type Stack + integer(IntKi) :: level + logical :: processed = .false. + real(DbKi) :: pos1(3), posMid(3), pos2(3) + real(DbKi) :: r1, rMid, r2 + real(DbKi) :: dl + integer(IntKi) :: secStat1, secStatMid, secStat2 + real(DbKi) :: dFdl1(6), dFdlMid(6), dFdl2(6) + end type + + type (Stack) :: SA(2*maxRecurLvl) ! Stack array + type (Stack) :: SE ! Stack element + ErrStat = ErrID_None - posMidL = 0.5_DbKi*(pos1+posMid) - posMidR = 0.5_DbKi*(posMid+pos2) - rMidL = 0.5_DbKi*(r1+rMid) - rMidR = 0.5_DbKi*(rMid+r2) + ! Convert unit vectors to ReKi for compatibility with GetSectionFreeSurfaceIntersects_Cyl + k_hat_Re = real(k_hat, ReKi) + y_hat_Re = real(y_hat, ReKi) + z_hat_Re = real(z_hat, ReKi) + n_hat_Re = real(n_hat, ReKi) + + ! Initialize first stack element + SA(1)%level = 1 + SA(1)%processed = .false. + SA(1)%pos1 = pos1 + SA(1)%posMid = posMid + SA(1)%pos2 = pos2 + SA(1)%r1 = r1 + SA(1)%rMid = rMid + SA(1)%r2 = r2 + SA(1)%dFdl1 = dFdl1 + SA(1)%dFdlMid = dFdlMid + SA(1)%dFdl2 = dFdl2 + SA(1)%dl = dl + SA(1)%secStat1 = secStat1 + SA(1)%secStatMid = secStatMid + SA(1)%secStat2 = secStat2 + + ! Initalize stack index + i = 1 + + ! Initialize output force + F_B_5pt = 0.0_DbKi + + ! Loop until all stack elements have been processed + do while (.true.) + + ! Find the next unprocessed stack element or exit if all processed + do while (SA(i)%processed) + i = i - 1 + if (i == 0) return ! All stack elements processed + end do - ! Avoid sections coincident with the SWL - IF ( ABS(k_hat(3)) > 0.999999_ReKi ) THEN ! Vertical member - IF ( EqualRealNos( posMidL(3), 0.0_DbKi ) ) THEN - posMidL(3) = posMidL(3) - 1.0E-6 * dl - END IF - IF ( EqualRealNos( posMidR(3), 0.0_DbKi ) ) THEN - posMidR(3) = posMidR(3) - 1.0E-6 * dl - END IF - END IF + ! Copy current stack element to local variable for processing + SE = SA(i) - ! Total hydrostatic load on the element (Simpsons Rule) - F_B_3pt = (dFdl1 + 4.0_DbKi*dFdlMid + dFdl2) * dl/6.0_DbKi + ! Compute mid points of left and right sub-elements + posMidL = 0.5_DbKi*(SE%pos1 + SE%posMid) + posMidR = 0.5_DbKi*(SE%posMid + SE%pos2) + rMidL = 0.5_DbKi*(SE%r1 + SE%rMid) + rMidR = 0.5_DbKi*(SE%rMid + SE%r2) - ! Mid point of left section - CALL GetSectionFreeSurfaceIntersects_Cyl( posMidL, FSPt, REAL(k_hat,ReKi), REAL(y_hat,ReKi), REAL(z_hat,ReKi), REAL(n_hat,ReKi), rMidL, theta1, theta2, secStatMidL) - CALL GetSectionHstLds_Cyl(p, origin, posMidL, k_hat, y_hat, z_hat, rMidL, dRdl, theta1, theta2, dFdlMidL) + ! Avoid sections coincident with the SWL + IF ( ABS(k_hat(3)) > 0.999999_ReKi ) THEN ! Vertical member + IF ( EqualRealNos( posMidL(3), 0.0_DbKi ) ) THEN + posMidL(3) = posMidL(3) - 1.0E-6 * SE%dl + END IF + IF ( EqualRealNos( posMidR(3), 0.0_DbKi ) ) THEN + posMidR(3) = posMidR(3) - 1.0E-6 * SE%dl + END IF + END IF - ! Mid point of right section - CALL GetSectionFreeSurfaceIntersects_Cyl( posMidR, FSPt, REAL(k_hat,ReKi), REAL(y_hat,ReKi), REAL(z_hat,ReKi), REAL(n_hat,ReKi), rMidR, theta1, theta2, secStatMidR) - CALL GetSectionHstLds_Cyl(p, origin, posMidR, k_hat, y_hat, z_hat, rMidR, dRdl, theta1, theta2, dFdlMidR) - - F_B_5pt = (dFdl1 + 4.0_DbKi*dFdlMidL + 2.0_DbKi*dFdlMid + 4.0_DbKi*dFdlMidR + dFdl2) * dl/12.0_DbKi + ! Total hydrostatic load on the element (3 point Simpsons Rule) + F_B_3pt_sub = (SE%dFdl1 + 4.0_DbKi*SE%dFdlMid + SE%dFdl2) * SE%dl/6.0_DbKi - ! Calculate error and check against tolerance - error = ABS(F_B_3pt - F_B_5pt) - tolMet = all(error <= MAX(RelTol*ABS(F_B_5pt),AbsTol)) + ! Mid point of left section + CALL GetSectionFreeSurfaceIntersects_Cyl( posMidL, FSPt, k_hat_Re, y_hat_Re, z_hat_Re, n_hat_Re, rMidL, theta1, theta2, secStatMidL) + CALL GetSectionHstLds_Cyl(p, origin, posMidL, k_hat, y_hat, z_hat, rMidL, dRdl, theta1, theta2, dFdlMidL) - ! If tolerance was met and (sub)element does not bound the waterplane, exit recursion - if (tolMet .and. (ABS(secStat1-secStat2) <= 1)) return + ! Mid point of right section + CALL GetSectionFreeSurfaceIntersects_Cyl( posMidR, FSPt, k_hat_Re, y_hat_Re, z_hat_Re, n_hat_Re, rMidR, theta1, theta2, secStatMidR) + CALL GetSectionHstLds_Cyl(p, origin, posMidR, k_hat, y_hat, z_hat, rMidR, dRdl, theta1, theta2, dFdlMidR) + + ! Total hydrostatic load on the element (5 point Simpsons Rule) + F_B_5pt_sub = (SE%dFdl1 + 4.0_DbKi*dFdlMidL + 2.0_DbKi*SE%dFdlMid + 4.0_DbKi*dFdlMidR + SE%dFdl2) * SE%dl/12.0_DbKi - ! If maximum recursion level reached, exit recursion with warning - if (recurLvl > maxRecurLvl) then - ErrStat = ErrID_Warn - return - end if + ! Calculate error and check against tolerance + error = ABS(F_B_3pt_sub - F_B_5pt_sub) + tolMet = all(error <= MAX(RelTol*ABS(F_B_5pt_sub),AbsTol)) + + ! If tolerance was met and (sub)element does not bound the waterplane, + ! Set processed flag, sum force, and continue + if (tolMet .and. (ABS(SE%secStat1 - SE%secStat2) <= 1)) then + SA(i)%processed = .true. + F_B_5pt = F_B_5pt + F_B_5pt_sub + cycle + end if + + ! If recursinon limit reached, + ! Set procesed flag, set error flag, and continue + if (SE%level + 1 > maxRecurLvl) then + SA(i)%processed = .true. + ErrStat = ErrID_Warn + cycle + end if - ! Recursively refine the load integration if tolerance not met - CALL RefineElementHstLds_Cyl(p,origin,pos1,posMidL,posMid,FSPt,r1,rMidL,rMid,0.5_DbKi*dl,dRdl,secStat1,secStatMidL,secStatMid,k_hat,y_hat,z_hat,n_hat,dFdl1,dFdlMidL,dFdlMid, recurLvl+1, tmp, ErrStat) - CALL RefineElementHstLds_Cyl(p,origin,posMid,posMidR,pos2,FSPt,rMid,rMidR,r2,0.5_DbKi*dl,dRdl,secStatMid,secStatMidR,secStat2,k_hat,y_hat,z_hat,n_hat,dFdlMid,dFdlMidR,dFdl2, recurLvl+1, F_B_5pt, ErrStat) - F_B_5pt = F_B_5pt + tmp + ! Push new branches onto stack + SA(i)%level = SE%level + 1 + SA(i)%processed = .false. + SA(i)%pos1 = SE%pos1 + SA(i)%posMid = posMidL + SA(i)%pos2 = SE%posMid + SA(i)%r1 = SE%r1 + SA(i)%rMid = rMidL + SA(i)%r2 = SE%rMid + SA(i)%dl = 0.5_DbKi * SE%dl + SA(i)%secStat1 = SE%secStat1 + SA(i)%secStatMid = secStatMidL + SA(i)%secStat2 = SE%secStatMid + SA(i)%dFdl1 = SE%dFdl1 + SA(i)%dFdlMid = dFdlMidL + SA(i)%dFdl2 = SE%dFdlMid + + SA(i+1)%level = SE%level + 1 + SA(i+1)%processed = .false. + SA(i+1)%pos1 = SE%posMid + SA(i+1)%posMid = posMidR + SA(i+1)%pos2 = SE%pos2 + SA(i+1)%r1 = SE%rMid + SA(i+1)%rMid = rMidR + SA(i+1)%r2 = SE%r2 + SA(i+1)%dl = 0.5_DbKi * SE%dl + SA(i+1)%secStat1 = SE%secStatMid + SA(i+1)%secStatMid = secStatMidR + SA(i+1)%secStat2 = SE%secStat2 + SA(i+1)%dFdl1 = SE%dFdlMid + SA(i+1)%dFdlMid = dFdlMidR + SA(i+1)%dFdl2 = SE%dFdl2 + + ! Increment stack index + i = i + 1 + end do END SUBROUTINE RefineElementHstLds_Cyl From 9715b1ea76772ceff2a85470c500dc71ce46de91 Mon Sep 17 00:00:00 2001 From: Derek Slaughter Date: Wed, 31 Dec 2025 14:50:31 -0500 Subject: [PATCH 09/12] Minor performance improvements in OLAF --- modules/aerodyn/src/FVW_BiotSavart.f90 | 92 ++++++++++++++------------ modules/aerodyn/src/FVW_Wings.f90 | 16 +++-- 2 files changed, 58 insertions(+), 50 deletions(-) diff --git a/modules/aerodyn/src/FVW_BiotSavart.f90 b/modules/aerodyn/src/FVW_BiotSavart.f90 index 038a38a953..fcc301511d 100644 --- a/modules/aerodyn/src/FVW_BiotSavart.f90 +++ b/modules/aerodyn/src/FVW_BiotSavart.f90 @@ -50,55 +50,61 @@ subroutine ui_seg_11(DeltaPa, DeltaPb, SegGamma, RegFunction, RegParam1, Uind) real(ReKi) :: xa, ya, za, xb, yb, zb !< Coordinates of X-Xa and X-Xb real(ReKi) :: exp_value !< ! - Uind(1:3)=0.0_ReKi + Uind(1:3) = 0.0_ReKi xa=DeltaPa(1); ya=DeltaPa(2); za=DeltaPa(3) xb=DeltaPb(1); yb=DeltaPb(2); zb=DeltaPb(3) norm_a = sqrt(xa*xa + ya*ya + za*za) norm_b = sqrt(xb*xb + yb*yb + zb*zb) denominator = norm_a*norm_b*(norm_a*norm_b + xa*xb+ya*yb+za*zb) ! |r1|*|r2|*(|r1|*|r2| + r1.r2) - if (denominator>PRECISION_UI) then - crossprod(1) = ya*zb-za*yb; crossprod(2) = za*xb-xa*zb; crossprod(3) = xa*yb-ya*xb - norm2_orth = crossprod(1)**2 + crossprod(2)**2 + crossprod(3)**2 - if (norm2_orth>PRECISION_UI) then ! On the singularity, Uind(1:3)=0.0_ReKi - norm2_r0 = (xa-xb)*(xa-xb) + (ya-yb)*(ya-yb) +(za-zb)*(za-zb) - if (norm2_r0>PRECISION_UI) then ! segment of zero length - ! --- Far field TODO - ! --- Regularization (close field) - norm2_orth = norm2_orth/norm2_r0 ! d = (r1xr2)/r0 - select case (RegFunction) ! - case ( idRegNone ) ! No vortex core model - Kv=1.0_ReKi - case ( idRegRankine ) ! Rankine - r_bar2 = norm2_orth/ RegParam1**2 - if (r_bar2<1) then - Kv=r_bar2 - else - Kv=1.0_ReKi - end if - case ( idRegLambOseen ) ! Lamb-Oseen - r_bar2 = norm2_orth/ RegParam1**2 - exp_value = -1.25643_ReKi*r_bar2 - if(exp_value Date: Wed, 31 Dec 2025 15:40:23 -0500 Subject: [PATCH 10/12] Apply suggestions from code review Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- modules/hydrodyn/src/Morison.f90 | 6 +++--- modules/moordyn/src/MoorDyn_Misc.f90 | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/modules/hydrodyn/src/Morison.f90 b/modules/hydrodyn/src/Morison.f90 index a28cc215e7..9baef40c9a 100644 --- a/modules/hydrodyn/src/Morison.f90 +++ b/modules/hydrodyn/src/Morison.f90 @@ -5233,7 +5233,7 @@ SUBROUTINE RefineElementHstLds_Cyl(p, origin, pos1, posMid, pos2, FSPt, r1, rMid REAL(DbKi) :: F_B_3pt_sub(6), F_B_5pt_sub(6) REAL(DbKi) :: error(6), tmp(6) LOGICAL :: tolMet - INTEGER(IntKi) :: i, ib + INTEGER(IntKi) :: i INTEGER(IntKi) :: secStatMidL, secStatMidR REAL(ReKi) :: k_hat_Re(3) REAL(ReKi) :: y_hat_Re(3) @@ -5342,8 +5342,8 @@ SUBROUTINE RefineElementHstLds_Cyl(p, origin, pos1, posMid, pos2, FSPt, r1, rMid cycle end if - ! If recursinon limit reached, - ! Set procesed flag, set error flag, and continue + ! If recursion limit reached, + ! Set processed flag, set error flag, and continue if (SE%level + 1 > maxRecurLvl) then SA(i)%processed = .true. ErrStat = ErrID_Warn diff --git a/modules/moordyn/src/MoorDyn_Misc.f90 b/modules/moordyn/src/MoorDyn_Misc.f90 index 03aed01572..9ae7a4fb6b 100644 --- a/modules/moordyn/src/MoorDyn_Misc.f90 +++ b/modules/moordyn/src/MoorDyn_Misc.f90 @@ -1010,7 +1010,7 @@ SUBROUTINE getWaterKin(p, m, x, y, z, t, U, Ud, zeta, PDyn, ErrStat, ErrMsg) ! SeaState wave kinematics case (3) - ! disable wavekin 3 during IC_gen, otherwise will never find steady state (becasue of waves) + ! disable wavekin 3 during IC_gen, otherwise will never find steady state (because of waves) if (m%IC_gen) return ! SeaState throws warning when queried location is out of bounds from the SeaState grid, so no need to handle here @@ -1019,7 +1019,7 @@ SUBROUTINE getWaterKin(p, m, x, y, z, t, U, Ud, zeta, PDyn, ErrStat, ErrMsg) ! (only pos needed because time is double in wave field, all other are outputs that will be set by WaveField_GetNodeWaveKin) xyz_sp = REAL([x, y, z], SiKi) - ! for now we will force the node to be in the water (forceNodeInWater = True). Rods handle partial submergence seperately so they need to get information from SeaState + ! for now we will force the node to be in the water (forceNodeInWater = True). Rods handle partial submergence separately so they need to get information from SeaState CALL WaveField_GetNodeWaveKin(p%WaveField, m%WaveField_m, t, xyz_sp, .true., .true., nodeInWater, WaveElev1, WaveElev2, zeta_sp, PDyn_sp, U_sp, Ud_sp, FAMCF, ErrStat2, ErrMsg2 ) ! outputs: nodeInWater, WaveElev1, WaveElev2, FAMCF all unused CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) From d14cdf4077a4db9e993ff30b25e1799dd715e3b7 Mon Sep 17 00:00:00 2001 From: Derek Slaughter Date: Wed, 31 Dec 2025 15:46:02 -0500 Subject: [PATCH 11/12] Add stack size check in RefineElementHstLds_Cyl --- modules/hydrodyn/src/Morison.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/modules/hydrodyn/src/Morison.f90 b/modules/hydrodyn/src/Morison.f90 index 9baef40c9a..f7f1153241 100644 --- a/modules/hydrodyn/src/Morison.f90 +++ b/modules/hydrodyn/src/Morison.f90 @@ -5342,9 +5342,9 @@ SUBROUTINE RefineElementHstLds_Cyl(p, origin, pos1, posMid, pos2, FSPt, r1, rMid cycle end if - ! If recursion limit reached, + ! If recursion limit reached or stack full, ! Set processed flag, set error flag, and continue - if (SE%level + 1 > maxRecurLvl) then + if ((SE%level + 1 > maxRecurLvl) .or. (i + 1 > size(SA))) then SA(i)%processed = .true. ErrStat = ErrID_Warn cycle From 11b4878e6b20d95f719e8104d40f05b489993cd1 Mon Sep 17 00:00:00 2001 From: Derek Slaughter Date: Mon, 5 Jan 2026 20:46:19 +0000 Subject: [PATCH 12/12] Remove recursion from RefineElementHstLds_Rec --- modules/hydrodyn/src/Morison.f90 | 195 ++++++++++++++++++++++--------- 1 file changed, 142 insertions(+), 53 deletions(-) diff --git a/modules/hydrodyn/src/Morison.f90 b/modules/hydrodyn/src/Morison.f90 index f7f1153241..096f3767cf 100644 --- a/modules/hydrodyn/src/Morison.f90 +++ b/modules/hydrodyn/src/Morison.f90 @@ -5389,7 +5389,7 @@ SUBROUTINE RefineElementHstLds_Cyl(p, origin, pos1, posMid, pos2, FSPt, r1, rMid END SUBROUTINE RefineElementHstLds_Cyl - RECURSIVE SUBROUTINE RefineElementHstLds_Rec(p, origin, pos1, posMid, pos2, FSPt, Sa1, SaMid, Sa2, Sb1, SbMid, Sb2, dl, dSadl, dSbdl, & + SUBROUTINE RefineElementHstLds_Rec(p, origin, pos1, posMid, pos2, FSPt, Sa1, SaMid, Sa2, Sb1, SbMid, Sb2, dl, dSadl, dSbdl, & secStat1, secStatMid, secStat2, k_hat, x_hat, y_hat, n_hat, dFdl1, dFdlMid, dFdl2, recurLvl, F_B_5pt, ErrStat, ErrMsg) TYPE(Morison_ParameterType), INTENT( IN ) :: p REAL(DbKi), INTENT( IN ) :: origin(3) @@ -5419,9 +5419,10 @@ RECURSIVE SUBROUTINE RefineElementHstLds_Rec(p, origin, pos1, posMid, pos2, FSPt REAL(DbKi) :: posMidL(3), posMidR(3) REAL(DbKi) :: SaMidL, SbMidL, SaMidR, SbMidR - REAL(DbKi) :: dFdlMidL(6), dFdlMidR(6), F_B_3pt(6) + REAL(DbKi) :: dFdlMidL(6), dFdlMidR(6) + REAL(DbKi) :: F_B_3pt_sub(6), F_B_5pt_sub(6) REAL(DbKi) :: error(6), tmp(6) - LOGICAL :: refine, tolMet + LOGICAL :: tolMet INTEGER(IntKi) :: i INTEGER(IntKi) :: secStatMidL, secStatMidR REAL(DbKi), PARAMETER :: RelTol = 1.0E-6 @@ -5429,64 +5430,152 @@ RECURSIVE SUBROUTINE RefineElementHstLds_Rec(p, origin, pos1, posMid, pos2, FSPt INTEGER(IntKi), PARAMETER :: maxRecurLvl = 50 CHARACTER(*), PARAMETER :: RoutineName = "RefineElementHstLds_Rec" + type Stack + integer(IntKi) :: level + logical :: processed = .false. + real(DbKi) :: pos1(3), posMid(3), pos2(3) + real(DbKi) :: Sa1, SaMid, Sa2 + real(DbKi) :: Sb1, SbMid, Sb2 + real(DbKi) :: dl + integer(IntKi) :: secStat1, secStatMid, secStat2 + real(DbKi) :: dFdl1(6), dFdlMid(6), dFdl2(6) + end type + + type (Stack) :: SA(2*maxRecurLvl) ! Stack array + type (Stack) :: SE ! Stack element + ErrStat = ErrID_None - ErrMsg = "" - posMidL = 0.5_DbKi*(pos1+posMid) - posMidR = 0.5_DbKi*(posMid+pos2) - SaMidL = 0.5_DbKi*(Sa1+SaMid) - SbMidL = 0.5_DbKi*(Sb1+SbMid) - SaMidR = 0.5_DbKi*(SaMid+Sa2) - SbMidR = 0.5_DbKi*(SbMid+Sb2) + ! Initialize first stack element + SA(1)%level = 1 + SA(1)%processed = .false. + SA(1)%pos1 = pos1 + SA(1)%posMid = posMid + SA(1)%pos2 = pos2 + SA(1)%Sa1 = Sa1 + SA(1)%SaMid = SaMid + SA(1)%Sa2 = Sa2 + SA(1)%Sb1 = Sb1 + SA(1)%SbMid = SbMid + SA(1)%Sb2 = Sb2 + SA(1)%dFdl1 = dFdl1 + SA(1)%dFdlMid = dFdlMid + SA(1)%dFdl2 = dFdl2 + SA(1)%dl = dl + SA(1)%secStat1 = secStat1 + SA(1)%secStatMid = secStatMid + SA(1)%secStat2 = secStat2 - ! Avoid sections coincident with the SWL - IF ( ABS(k_hat(3)) > 0.999999_ReKi ) THEN ! Vertical member - IF ( EqualRealNos( posMidL(3), 0.0_DbKi ) ) THEN - posMidL(3) = posMidL(3) - 1.0E-6 * dl - END IF - IF ( EqualRealNos( posMidR(3), 0.0_DbKi ) ) THEN - posMidR(3) = posMidR(3) - 1.0E-6 * dl - END IF - END IF + ! Initalize stack index + i = 1 - ! Total hydrostatic load on the element (Simpsons Rule) - F_B_3pt = (dFdl1 + 4.0_DbKi*dFdlMid + dFdl2) * dl/6.0_DbKi + ! Initialize output force + F_B_5pt = 0.0_DbKi - ! Mid point of left section - CALL GetSectionHstLds_Rec(p, origin, posMidL, k_hat, x_hat, y_hat, SaMidL, SbMidL, dSadl, dSbdl, FSPt, n_hat, dFdlMidL, secStatMidL) + ! Loop until all stack elements have been processed + do while (.true.) - ! Mid point of right section - CALL GetSectionHstLds_Rec(p, origin, posMidR, k_hat, x_hat, y_hat, SaMidR, SbMidR, dSadl, dSbdl, FSPt, n_hat, dFdlMidR, secStatMidR) - - F_B_5pt = (dFdl1 + 4.0_DbKi*dFdlMidL + 2.0_DbKi*dFdlMid + 4.0_DbKi*dFdlMidR + dFdl2) * dl/12.0_DbKi + ! Find the next unprocessed stack element or exit if all processed + do while (SA(i)%processed) + i = i - 1 + if (i == 0) return ! All stack elements processed + end do - error = ABS(F_B_3pt - F_B_5pt) - tolMet = .TRUE. - DO i = 1,6 - IF ( error(i) > MAX(RelTol*ABS(F_B_5pt(i)),AbsTol) ) THEN - tolMet = .FALSE. - END IF - END DO - refine = .NOT. tolMet - IF (ABS(secStat1-secStat2)>1) THEN ! (Sub)element bounds the waterplane - refine = .TRUE. ! Keep refining irrespective of tolMet to avoid premature termination - END IF - IF ( recurLvl > maxRecurLvl ) THEN - refine = .FALSE. - IF (.NOT. tolMet) THEN - CALL SetErrStat(ErrID_Warn, 'Tolerance for element hydrostatic load not met after the maximum allowed level of recursion is reached. Consider reducing MDivSize.', ErrStat, ErrMsg, RoutineName ) - ! ELSE - ! Free surface is likely normal to the element. + ! Copy current stack element to local variable for processing + SE = SA(i) + + ! Compute mid points of left and right sub-elements + posMidL = 0.5_DbKi*(SE%pos1+SE%posMid) + posMidR = 0.5_DbKi*(SE%posMid+SE%pos2) + SaMidL = 0.5_DbKi*(SE%Sa1+SE%SaMid) + SbMidL = 0.5_DbKi*(SE%Sb1+SE%SbMid) + SaMidR = 0.5_DbKi*(SE%SaMid+SE%Sa2) + SbMidR = 0.5_DbKi*(SE%SbMid+SE%Sb2) + + ! Avoid sections coincident with the SWL + IF ( ABS(k_hat(3)) > 0.999999_ReKi ) THEN ! Vertical member + IF ( EqualRealNos( posMidL(3), 0.0_DbKi ) ) THEN + posMidL(3) = posMidL(3) - 1.0E-6 * SE%dl + END IF + IF ( EqualRealNos( posMidR(3), 0.0_DbKi ) ) THEN + posMidR(3) = posMidR(3) - 1.0E-6 * SE%dl + END IF END IF - END IF - - IF (refine) THEN ! Recursively refine the load integration if tolerance not met - CALL RefineElementHstLds_Rec(p,origin,pos1,posMidL,posMid,FSPt,Sa1,SaMidL,SaMid,Sb1,SbMidL,SbMid,0.5_DbKi*dl,dSadl,dSbdl, & - secStat1,secStatMidL,secStatMid,k_hat,x_hat,y_hat,n_hat,dFdl1,dFdlMidL,dFdlMid,recurLvl+1,tmp,ErrStat,ErrMsg) - CALL RefineElementHstLds_Rec(p,origin,posMid,posMidR,pos2,FSPt,SaMid,SaMidR,Sa2,SbMid,SbMidR,Sb2,0.5_DbKi*dl,dSadl,dSbdl, & - secStatMid,secStatMidR,secStat2,k_hat,x_hat,y_hat,n_hat,dFdlMid,dFdlMidR,dFdl2,recurLvl+1,F_B_5pt,ErrStat,ErrMsg) - F_B_5pt = F_B_5pt + tmp - END IF + + ! Total hydrostatic load on the element (Simpsons Rule) + F_B_3pt_sub = (SE%dFdl1 + 4.0_DbKi*SE%dFdlMid + SE%dFdl2) * SE%dl/6.0_DbKi + + ! Mid point of left section + CALL GetSectionHstLds_Rec(p, origin, posMidL, k_hat, x_hat, y_hat, SaMidL, SbMidL, dSadl, dSbdl, FSPt, n_hat, dFdlMidL, secStatMidL) + + ! Mid point of right section + CALL GetSectionHstLds_Rec(p, origin, posMidR, k_hat, x_hat, y_hat, SaMidR, SbMidR, dSadl, dSbdl, FSPt, n_hat, dFdlMidR, secStatMidR) + + ! Total hydrostatic load on the element (5 point Simpsons Rule) + F_B_5pt_sub = (SE%dFdl1 + 4.0_DbKi*dFdlMidL + 2.0_DbKi*SE%dFdlMid + 4.0_DbKi*dFdlMidR + SE%dFdl2) * SE%dl/12.0_DbKi + + ! Calculate error and check against tolerance + error = ABS(F_B_3pt_sub - F_B_5pt_sub) + tolMet = all(error <= MAX(RelTol*ABS(F_B_5pt_sub),AbsTol)) + + ! If tolerance was met and (sub)element does not bound the waterplane, + ! Set processed flag, sum force, and continue + if (tolMet .and. (ABS(SE%secStat1 - SE%secStat2) <= 1)) then + SA(i)%processed = .true. + F_B_5pt = F_B_5pt + F_B_5pt_sub + cycle + end if + + ! If recursion limit reached or stack full, + ! Set processed flag, set error flag, and continue + if ((SE%level + 1 > maxRecurLvl) .or. (i + 1 > size(SA))) then + SA(i)%processed = .true. + ErrStat = ErrID_Warn + cycle + end if + + ! Push new branches onto stack + SA(i)%level = SE%level + 1 + SA(i)%processed = .false. + SA(i)%pos1 = SE%pos1 + SA(i)%posMid = posMidL + SA(i)%pos2 = SE%posMid + SA(i)%Sa1 = SE%Sa1 + SA(i)%SaMid = SaMidL + SA(i)%Sa2 = SE%SaMid + SA(i)%Sb1 = SE%Sb1 + SA(i)%SbMid = SbMidL + SA(i)%Sb2 = SE%SbMid + SA(i)%dl = 0.5_DbKi * SE%dl + SA(i)%secStat1 = SE%secStat1 + SA(i)%secStatMid = secStatMidL + SA(i)%secStat2 = SE%secStatMid + SA(i)%dFdl1 = SE%dFdl1 + SA(i)%dFdlMid = dFdlMidL + SA(i)%dFdl2 = SE%dFdlMid + + SA(i+1)%level = SE%level + 1 + SA(i+1)%processed = .false. + SA(i+1)%pos1 = SE%posMid + SA(i+1)%posMid = posMidR + SA(i+1)%pos2 = SE%pos2 + SA(i+1)%Sa1 = SE%SaMid + SA(i+1)%SaMid = SaMidR + SA(i+1)%Sa2 = SE%Sa2 + SA(i+1)%Sb1 = SE%SbMid + SA(i+1)%SbMid = SbMidR + SA(i+1)%Sb2 = SE%Sb2 + SA(i+1)%dl = 0.5_DbKi * SE%dl + SA(i+1)%secStat1 = SE%secStatMid + SA(i+1)%secStatMid = secStatMidR + SA(i+1)%secStat2 = SE%secStat2 + SA(i+1)%dFdl1 = SE%dFdlMid + SA(i+1)%dFdlMid = dFdlMidR + SA(i+1)%dFdl2 = SE%dFdl2 + + ! Increment stack index + i = i + 1 + end do END SUBROUTINE RefineElementHstLds_Rec