diff --git a/modules/aerodyn/src/FVW_BiotSavart.f90 b/modules/aerodyn/src/FVW_BiotSavart.f90 index 36c73d1dd9..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_valuePRECISION_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 +321,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 diff --git a/modules/aerodyn/src/FVW_Wings.f90 b/modules/aerodyn/src/FVW_Wings.f90 index 443653caac..4279d9d59c 100644 --- a/modules/aerodyn/src/FVW_Wings.f90 +++ b/modules/aerodyn/src/FVW_Wings.f90 @@ -281,23 +281,25 @@ subroutine Wings_ComputeCirculation(t, z, z_prev, p, x, m, AFInfo, ErrStat, ErrM GammaScale=1.0_ReKi endif - if (p%CircSolvMethod==idCircPrescribed) then + select case (p%CircSolvMethod) + ! Prescribed circulation + case (idCircPrescribed) do iW = 1, p%nWings !Loop over lifting lines z%W(iW)%Gamma_LL(1:p%W(iW)%nSpan) = p%W(iW)%PrescribedCirculation(1:p%W(iW)%nSpan) m%W(iW)%Vind_CP=-9999._ReKi !< Safety m%W(iW)%Vtot_CP=-9999._ReKi !< Safety enddo - else if (p%CircSolvMethod==idCircPolarData) then - ! --- Solve for circulation using polar data + ! Solve for circulation using polar data + case (idCircPolarData) CALL Wings_ComputeCirculationPolarData(z, z_prev, p, x, m, AFInfo, GammaScale, ErrStat, ErrMsg, iLabel) - else if (p%CircSolvMethod==idCircNoFlowThrough) then + case (idCircNoFlowThrough) ! --- Solve for circulation using the no-flow through condition ErrMsg='Circulation method nor implemented'; ErrStat=ErrID_Fatal; return ! should never happen - else + case default ErrMsg='Circulation method nor implemented'; ErrStat=ErrID_Fatal; return ! should never happen - endif + end select ! Scale circulation (for initial transient) do iW = 1, p%nWings !Loop over lifting lines @@ -428,9 +430,9 @@ subroutine Wings_ComputeCirculationPolarData(z, z_prev, p, x, m, AFInfo, GammaSc do iDepth=1,p%iNWStart ! Two first panels ! --- Defining a ring P1=x%W(iW)%r_NW(1:3,iSpan ,iDepth ) + P4=x%W(iW)%r_NW(1:3,iSpan ,iDepth+1) P2=x%W(iW)%r_NW(1:3,iSpan+1,iDepth ) P3=x%W(iW)%r_NW(1:3,iSpan+1,iDepth+1) - P4=x%W(iW)%r_NW(1:3,iSpan ,iDepth+1) ! --- Induced velocity from ring, on all other control points (have to loop on rotors and wings and span again) kCP2=1 do iWCP=1,p%nWings 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..096f3767cf 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) - + 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,84 +5223,175 @@ 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) :: 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(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 - ErrMsg = "" - 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 + ! 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) + rMidL = 0.5_DbKi*(SE%r1 + SE%rMid) + rMidR = 0.5_DbKi*(SE%rMid + SE%r2) + + ! 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 - ! Total hydrostatic load on the element (Simpsons Rule) - F_B_3pt = (dFdl1 + 4.0_DbKi*dFdlMid + dFdl2) * dl/6.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 - ! 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) + ! 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) - ! 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) - - F_B_5pt = (dFdl1 + 4.0_DbKi*dFdlMidL + 2.0_DbKi*dFdlMid + 4.0_DbKi*dFdlMidR + dFdl2) * dl/12.0_DbKi + ! 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 - 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 + ! 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)%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 - RECURSIVE SUBROUTINE RefineElementHstLds_Rec( 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) REAL(DbKi), INTENT( IN ) :: pos1(3) REAL(DbKi), INTENT( IN ) :: posMid(3) @@ -5287,9 +5419,10 @@ RECURSIVE SUBROUTINE RefineElementHstLds_Rec( 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 @@ -5297,69 +5430,157 @@ RECURSIVE SUBROUTINE RefineElementHstLds_Rec( 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 + + ! 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 + ! 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 - ! Total hydrostatic load on the element (Simpsons Rule) - F_B_3pt = (dFdl1 + 4.0_DbKi*dFdlMid + dFdl2) * dl/6.0_DbKi + ! 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( origin, posMidL, k_hat, x_hat, y_hat, SaMidL, SbMidL, dSadl, dSbdl, FSPt, n_hat, dFdlMidL, secStatMidL) + ! 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( 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 + ! 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) - 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_Rec(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, & - 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 (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 - END SUBROUTINE RefineElementHstLds_Rec + ! 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 - SUBROUTINE GetEndPlateHstLds_Cyl(pos0, k_hat, y_hat, z_hat, R, theta1, theta2, F) + ! 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 + + 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 +5624,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 +5646,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 +5719,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 +5742,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 +5763,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 +5801,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 +5862,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 +5892,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 +5904,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 +5915,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 +5944,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 +6005,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 +6031,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 +6061,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 +6171,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 +6344,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) 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/moordyn/src/MoorDyn_Misc.f90 b/modules/moordyn/src/MoorDyn_Misc.f90 index 65fd5edddd..9ae7a4fb6b 100644 --- a/modules/moordyn/src/MoorDyn_Misc.f90 +++ b/modules/moordyn/src/MoorDyn_Misc.f90 @@ -948,28 +948,21 @@ 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 outputs 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) + ! no wave kinematics, do nothing (already initialized to zero above) + case (0) - ! 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) - - ELSEIF (p%WaterKin == 1 .OR. p%WaterKin == 2) THEN ! old or hybrid approach. SeaState contributions handeled in setupWaterKin, just proceed using old method + ! old or hybrid approach. SeaState contributions handeled in setupWaterKin, just proceed using old method + case (1,2) ! If wave kinematics enabled, get interpolated values from grid - IF (p%WaveKin > 0) THEN + 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 @@ -977,59 +970,68 @@ SUBROUTINE getWaterKin(p, m, x, y, z, t, U, Ud, zeta, PDyn, ErrStat, ErrMsg) 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 + 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 ) + zp = (z - zeta) * p%WtrDpth/(p%WtrDpth + zeta) - CALL getInterpNumsSiKi(p%pzWave , REAL(zp,SiKi), 1, iz, fz) ! wave grid + 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) ) - - 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 - - ENDIF + 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 - 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 + ! SeaState wave kinematics + case (3) + + ! 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 + + ! Pack all MD inputs to WaveGrid input data types (double to single) + ! (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 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) + + ! 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) + + case default + call SetErrStat(ErrID_Fatal, "Invalid value of p%WaterKin", ErrStat, ErrMsg, RoutineName) + end select END SUBROUTINE getWaterKin 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_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 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/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) 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: