diff --git a/modules/moordyn/src/MoorDyn.f90 b/modules/moordyn/src/MoorDyn.f90 index ea63f639c6..ebde8cf519 100644 --- a/modules/moordyn/src/MoorDyn.f90 +++ b/modules/moordyn/src/MoorDyn.f90 @@ -200,14 +200,15 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er m%PtfmInit = InitInp%PtfmInit(:,1) ! is this copying necssary in case this is an individual instance in FAST.Farm? - - ! Check if this MoorDyn instance is being run from FAST.Farm (indicated by FarmSize > 0) if (InitInp%FarmSize > 0) then CALL WrScr(' >>> MoorDyn is running in array mode <<< ') ! could make sure the size of this is right: SIZE(InitInp%FarmCoupledKinematics) p%nTurbines = InitInp%FarmSize - else ! FarmSize==0 indicates normal, FAST module mode + else if (InitInp%FarmSize < 0) then ! Farmsize==-1 indicates standlone, run MoorDyn as a standalone code with no openfast coupling + p%Standalone = 1 + p%nTurbines = 1 + else ! FarmSize==0 indicates normal, FAST module mode p%nTurbines = 1 ! if a regular FAST module mode, we treat it like a nTurbine=1 farm case END IF @@ -244,29 +245,29 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! ----------------------------------------------------------------- ! Read the primary MoorDyn input file, or copy from passed input - if (InitInp%UsePrimaryInputFile) then - ! Read the entire input file, minus any comment lines, into the FileInfo_In - ! data structure in memory for further processing. - call ProcessComFile( InitInp%FileName, FileInfo_In, ErrStat2, ErrMsg2 ) - CALL GetPath( InitInp%FileName, p%PriPath ) ! Input files will be relative to the path where the primary input file is located. - else - call NWTC_Library_CopyFileInfoType( InitInp%PassedPrimaryInputData, FileInfo_In, MESH_NEWCOPY, ErrStat2, ErrMsg2 ) - p%PriPath = "" - endif - if (Failed()) return; - - ! For diagnostic purposes, the following can be used to display the contents - ! of the FileInfo_In data structure. - !call Print_FileInfo_Struct( CU, FileInfo_In ) ! CU is the screen -- different number on different systems. + if (InitInp%UsePrimaryInputFile) then + ! Read the entire input file, minus any comment lines, into the FileInfo_In + ! data structure in memory for further processing. + call ProcessComFile( InitInp%FileName, FileInfo_In, ErrStat2, ErrMsg2 ) + CALL GetPath( InitInp%FileName, p%PriPath ) ! Input files will be relative to the path where the primary input file is located. + else + call NWTC_Library_CopyFileInfoType( InitInp%PassedPrimaryInputData, FileInfo_In, MESH_NEWCOPY, ErrStat2, ErrMsg2 ) + p%PriPath = "" + endif + if (Failed()) return; + + ! For diagnostic purposes, the following can be used to display the contents + ! of the FileInfo_In data structure. + !call Print_FileInfo_Struct( CU, FileInfo_In ) ! CU is the screen -- different number on different systems. ! Parse the FileInfo_In structure of data from the inputfile into the InitInp%InputFile structure -! CALL ParsePrimaryFileInfo_BuildModel( PriPath, InitInp, FileInfo_In, InputFileDat, p, m, UnEc, ErrStat2, ErrMsg2 ) -! if (Failed()) return; + ! CALL ParsePrimaryFileInfo_BuildModel( PriPath, InitInp, FileInfo_In, InputFileDat, p, m, UnEc, ErrStat2, ErrMsg2 ) + ! if (Failed()) return; -!NOTE: This could be split into a separate routine for easier to read code + !NOTE: This could be split into a separate routine for easier to read code !------------------------------------------------------------------------------------------------- ! Parsing of input file from the FileInfo_In data structure ! - FileInfo_Type is essentially a string array with some metadata. @@ -424,6 +425,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er END IF write(p%UnLog,'(A)', IOSTAT=ErrStat2) "MoorDyn v2 log file with output level "//TRIM(Num2LStr(p%writeLog)) write(p%UnLog,'(A)', IOSTAT=ErrStat2) "Note: options above the writeLog line in the input file will not be recorded." + write(p%UnLog,'(A)', IOSTAT=ErrStat2) " Input File Summary:" end if else if ( OptString == 'DTM') THEN read (OptValue,*) p%dtM0 @@ -457,14 +459,36 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er read (OptValue,*) p%mc else if ( OptString == 'CV') then read (OptValue,*) p%cv + else if ( OptString == 'INERTIALF') then + read (OptValue,*) p%inertialF else CALL SetErrStat( ErrID_Warn, 'Unable to interpret input '//trim(OptString)//' in OPTIONS section.', ErrStat, ErrMsg, RoutineName ) end if nOpts = nOpts + 1 Line = NextLine(i) + END DO + if (p%writeLog > 1) then + write(p%UnLog, '(A)' ) " - Options List:" + write(p%UnLog, '(A17,f12.4)') " dtm : ", p%dtM0 + write(p%UnLog, '(A17,f12.4)') " g : ", p%g + write(p%UnLog, '(A17,f12.4)') " rhoW : ", p%rhoW + write(p%UnLog, '(A17,A)' ) " Depth : ", DepthValue ! water depth input read in as a string to be processed by setupBathymetry + write(p%UnLog, '(A17,f12.4)') " kBot : ", p%kBot + write(p%UnLog, '(A17,f12.4)') " cBot : ", p%cBot + write(p%UnLog, '(A17,f12.4)') " dtIC : ", InputFileDat%dtIC + write(p%UnLog, '(A17,f12.4)') " TMaxIC : ", InputFileDat%TMaxIC + write(p%UnLog, '(A17,f12.4)') " CdScaleIC: ", InputFileDat%CdScaleIC + write(p%UnLog, '(A17,f12.4)') " threshIC : ", InputFileDat%threshIC + write(p%UnLog, '(A17,A)' ) " WaterKin : ", WaterKinValue + write(p%UnLog, '(A17,f12.4)') " dtOut : ", p%dtOut + write(p%UnLog, '(A17,f12.4)') " mu_kT : ", p%mu_kT + write(p%UnLog, '(A17,f12.4)') " mu_kA : ", p%mu_kA + write(p%UnLog, '(A17,f12.4)') " mc : ", p%mc + write(p%UnLog, '(A17,f12.4)') " cv : ", p%cv + end if else if (INDEX(Line, "OUTPUT") > 0) then ! if output header @@ -583,8 +607,8 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er Line = NextLine(i) Line = NextLine(i) - ! process each line - DO l = 1,p%nLineTypes + ! process each line + DO l = 1,p%nLineTypes !read into a line Line = NextLine(i) @@ -636,6 +660,9 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er read(tempStrings(2), *) m%LineTypeList(l)%BA_D else if (m%LineTypeList(l)%ElasticMod == 2) then ! case where there is no dynamic damping for viscoelastic model (will it work)? CALL WrScr("Warning, viscoelastic model being used with zero damping on the dynamic stiffness.") + if (p%writeLog > 0) then + write(p%UnLog,'(A)') "Warning, viscoelastic model being used with zero damping on the dynamic stiffness." + end if end if ! get the regular/static coefficient or relation in all cases (can be from a lookup table?) CALL getCoefficientOrCurve(tempStrings(1), m%LineTypeList(l)%BA, & @@ -653,16 +680,16 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er m%LineTypeList(l)%IdNum = l ! write lineType information to log file - if (p%writeLog > 1) then - write(p%UnLog, '(A12,A20)' ) " LineType"//trim(num2lstr(l))//":" - write(p%UnLog, '(A12,A20)' ) " name: ", m%LineTypeList(l)%name - write(p%UnLog, '(A12,f12.4)') " d : ", m%LineTypeList(l)%d - write(p%UnLog, '(A12,f12.4)') " w : ", m%LineTypeList(l)%w - write(p%UnLog, '(A12,f12.4)') " Cdn : ", m%LineTypeList(l)%Cdn - write(p%UnLog, '(A12,f12.4)') " Can : ", m%LineTypeList(l)%Can - write(p%UnLog, '(A12,f12.4)') " Cdt : ", m%LineTypeList(l)%Cdt - write(p%UnLog, '(A12,f12.4)') " Cat : ", m%LineTypeList(l)%Cat - end if + if (p%writeLog > 1) then + write(p%UnLog, '(A)' ) " - LineType"//trim(num2lstr(l))//":" + write(p%UnLog, '(A12,A)' ) " name: ", trim(m%LineTypeList(l)%name) + write(p%UnLog, '(A12,f12.4)') " d : ", m%LineTypeList(l)%d + write(p%UnLog, '(A12,f12.4)') " w : ", m%LineTypeList(l)%w + write(p%UnLog, '(A12,f12.4)') " Cdn : ", m%LineTypeList(l)%Cdn + write(p%UnLog, '(A12,f12.4)') " Can : ", m%LineTypeList(l)%Can + write(p%UnLog, '(A12,f12.4)') " Cdt : ", m%LineTypeList(l)%Cdt + write(p%UnLog, '(A12,f12.4)') " Cat : ", m%LineTypeList(l)%Cat + end if IF ( ErrStat2 /= ErrID_None ) THEN CALL SetErrStat( ErrID_Fatal, ErrMsg2, ErrStat, ErrMsg, RoutineName ) @@ -707,16 +734,16 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! specify IdNum of rod type for error checking m%RodTypeList(l)%IdNum = l - ! write lineType information to log file + ! write rodType information to log file if (p%writeLog > 1) then - write(p%UnLog, '(A12,A20)' ) " RodType"//trim(num2lstr(l))//":" - write(p%UnLog, '(A12,A20)' ) " name: ", m%RodTypeList(l)%name - write(p%UnLog, '(A12,f12.4)') " d : ", m%RodTypeList(l)%d - write(p%UnLog, '(A12,f12.4)') " w : ", m%RodTypeList(l)%w - write(p%UnLog, '(A12,f12.4)') " Cdn : ", m%RodTypeList(l)%Cdn - write(p%UnLog, '(A12,f12.4)') " Can : ", m%RodTypeList(l)%Can - write(p%UnLog, '(A12,f12.4)') " Cdt : ", m%RodTypeList(l)%CdEnd - write(p%UnLog, '(A12,f12.4)') " Cat : ", m%RodTypeList(l)%CaEnd + write(p%UnLog, '(A)' ) " - RodType"//trim(num2lstr(l))//":" + write(p%UnLog, '(A14,A)' ) " name: ", trim(m%RodTypeList(l)%name) + write(p%UnLog, '(A14,f12.4)') " d : ", m%RodTypeList(l)%d + write(p%UnLog, '(A14,f12.4)') " w : ", m%RodTypeList(l)%w + write(p%UnLog, '(A14,f12.4)') " Cdn : ", m%RodTypeList(l)%Cdn + write(p%UnLog, '(A14,f12.4)') " Can : ", m%RodTypeList(l)%Can + write(p%UnLog, '(A14,f12.4)') " CdEnd : ", m%RodTypeList(l)%CdEnd + write(p%UnLog, '(A14,f12.4)') " CaEnd : ", m%RodTypeList(l)%CaEnd end if IF ( ErrStat2 /= ErrID_None ) THEN @@ -768,7 +795,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er else CALL SetErrStat( ErrID_Fatal, 'Body '//trim(Num2LStr(l))//' CG entry (col 10) must have 1 or 3 numbers.' , ErrStat, ErrMsg, RoutineName ) end if - ! process mements of inertia + ! process moments of inertia CALL SplitByBars(tempString3, N, tempStrings) if (N == 1) then ! if only one entry, use it for all directions READ(tempString3, *) m%BodyList(l)%BodyI(1) @@ -813,6 +840,10 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er CALL WrScr(' Unable to parse Body '//trim(Num2LStr(l))//' on row '//trim(Num2LStr(i))//' in input file.') ! Specific screen output because errors likely CALL WrScr(' Ensure row has all 13 columns needed in MDv2 input file (13th Dec 2021).') CALL SetErrStat( ErrID_Fatal, 'Failed to read bodies.' , ErrStat, ErrMsg, RoutineName ) + if (p%writeLog > 0) then + write(p%UnLog,'(A)') ' Unable to parse Body '//trim(Num2LStr(l))//' on row '//trim(Num2LStr(i))//' in input file.' + write(p%UnLog,'(A)') ' Ensure row has all 13 columns needed in MDv2 input file (13th Dec 2021).' + end if CALL CleanUp() RETURN END IF @@ -834,6 +865,20 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er m%CpldBodyIs(p%nCpldBodies(1),1) = l ! body initial position due to coupling will be adjusted later + + else if ((let1 == "VESSELPINNED") .or. (let1 == "VESPIN") .or. (let1 == "COUPLEDPINNED") .or. (let1 == "CPLDPIN")) then ! if a pinned coupled body, add to list and add + m%BodyList(l)%typeNum = 2 + + p%nCpldBodies(1)=p%nCpldBodies(1)+1 ! add + p%nFreeBodies =p%nFreeBodies+1 ! add this pinned body to the free list because it is half free + + m%BodyStateIs1(p%nFreeBodies) = Nx+1 + m%BodyStateIsN(p%nFreeBodies) = Nx+6 + Nx = Nx + 6 ! add 6 state variables for each pinned body + + m%CpldBodyIs(p%nCpldBodies(1),1) = l + m%FreeBodyIs(p%nFreeBodies) = l + ! TODO: add option for body coupling to different turbines in FAST.Farm <<< @@ -874,6 +919,19 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er CALL CleanUp() RETURN END IF + + ! write body information to log file + if (p%writeLog > 1) then + write(p%UnLog, '(A)' ) " - Body"//trim(num2lstr(l))//":" + write(p%UnLog, '(A14,I2)' ) " id : ", m%BodyList(l)%IdNum + write(p%UnLog, '(A14,A)' ) " attach: ", trim(tempString1) + write(p%UnLog, '(A14,f12.4)') " v : ", m%BodyList(l)%bodyV + write(p%UnLog, '(A14,f12.4)') " m : ", m%BodyList(l)%bodyM + write(p%UnLog, '(A14,A)' ) " I : ", trim(num2lstr(m%BodyList(l)%BodyI(1)))//", "//trim(num2lstr(m%BodyList(l)%BodyI(2)))//", "//trim(num2lstr(m%BodyList(l)%BodyI(3))) + write(p%UnLog, '(A14,A)' ) " rCG : ", trim(num2lstr(m%BodyList(l)%rCG(1)))//", "//trim(num2lstr(m%BodyList(l)%rCG(2)))//", "//trim(num2lstr(m%BodyList(l)%rCG(3))) + write(p%UnLog, '(A14,A)' ) " CdA : ", trim(num2lstr(m%BodyList(l)%BodyCdA(1)))//", "//trim(num2lstr(m%BodyList(l)%BodyCdA(2)))//", "//trim(num2lstr(m%BodyList(l)%BodyCdA(3)))//", "//trim(num2lstr(m%BodyList(l)%BodyCdA(4)))//", "//trim(num2lstr(m%BodyList(l)%BodyCdA(5)))//", "//trim(num2lstr(m%BodyList(l)%BodyCdA(6))) + write(p%UnLog, '(A14,A)' ) " Ca : ", trim(num2lstr(m%BodyList(l)%BodyCa(1)))//", "//trim(num2lstr(m%BodyList(l)%BodyCa(2)))//", "//trim(num2lstr(m%BodyList(l)%BodyCa(3)))//", "//trim(num2lstr(m%BodyList(l)%BodyCa(4)))//", "//trim(num2lstr(m%BodyList(l)%BodyCa(5)))//", "//trim(num2lstr(m%BodyList(l)%BodyCa(6))) + end if IF (wordy > 1) print *, "Set up body ", l, " of type ", m%BodyList(l)%typeNum @@ -1005,10 +1063,10 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! TODO: add option for body coupling to different turbines in FAST.Farm <<< - else if ((let1 == "CONNECT") .or. (let1 == "CON") .or. (let1 == "FREE")) then + else if ((let1 == "ROD") .or. (let1 == "R") .or. (let1 == "FREE")) then m%RodList(l)%typeNum = 0 - p%nFreeRods=p%nFreeRods+1 ! add this pinned rod to the free list because it is half free + p%nFreeRods=p%nFreeRods+1 m%RodStateIs1(p%nFreeRods) = Nx+1 m%RodStateIsN(p%nFreeRods) = Nx+12 @@ -1050,6 +1108,14 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! specify IdNum of line for error checking m%RodList(l)%IdNum = l + if (p%writeLog > 1) then + write(p%UnLog, '(A)' ) " - Rod"//trim(num2lstr(m%RodList(l)%IdNum))//":" + write(p%UnLog, '(A15,I2)' ) " ID : ", m%RodList(l)%IdNum + write(p%UnLog, '(A15,A)' ) " Type : ", trim(m%RodTypeList(m%RodList(l)%PropsIdNum)%name) + write(p%UnLog, '(A15,A)' ) " Attach : ", trim(tempString2) + write(p%UnLog, '(A15,I2)' ) " NumSegs: ", m%RodList(l)%N + end if + ! check for sequential IdNums IF ( m%RodList(l)%IdNum .NE. l ) THEN CALL SetErrStat( ErrID_Fatal, 'Line numbers must be sequential starting from 1.', ErrStat, ErrMsg, RoutineName ) @@ -1105,6 +1171,9 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er if ((INDEX(tempString4, "SEABED") > 0 ) .or. (INDEX(tempString4, "GROUND") > 0 ) .or. (INDEX(tempString4, "FLOOR") > 0 )) then ! if keyword used CALL WrScr('Point '//trim(Num2LStr(l))//' depth set to be on the seabed; finding z location based on depth/bathymetry') ! interpret the anchor depth value as a 'seabed' input + if (p%writeLog > 0) then + write(p%UnLog,'(A)') 'Point '//trim(Num2LStr(l))//' depth set to be on the seabed; finding z location based on depth/bathymetry' + end if CALL getDepthFromBathymetry(m%BathymetryGrid, m%BathGrid_Xs, m%BathGrid_Ys, tempArray(1), tempArray(2), depth, nvec) ! meaning the anchor should be at the depth of the local bathymetry tempArray(3) = -depth else ! if the anchor depth input isn't one of the supported keywords, @@ -1124,6 +1193,10 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er CALL WrScr(' Unable to parse Point '//trim(Num2LStr(l))//' row in input file.') ! Specific screen output because errors likely CALL WrScr(' Ensure row has all 9 columns, including CdA and Ca.') ! to be caused by non-updated input file formats. CALL SetErrStat( ErrID_Fatal, 'Failed to read points.' , ErrStat, ErrMsg, RoutineName ) ! would be nice to specify which line <<<<<<<<< + if (p%writeLog > 0) then + write(p%UnLog,'(A)') ' Unable to parse Point '//trim(Num2LStr(l))//' row in input file.' + write(p%UnLog,'(A)') ' Ensure row has all 9 columns, including CdA and Ca.' + end if CALL CleanUp() RETURN END IF @@ -1188,6 +1261,9 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er p%nCpldPoints(J) = p%nCpldPoints(J) + 1 ! increment counter for the appropriate turbine m%CpldPointIs(p%nCpldPoints(J),J) = l CALL WrScr(' added point '//TRIM(int2lstr(l))//' as fairlead for turbine '//trim(int2lstr(J))) + if (p%writeLog > 0) then + write(p%UnLog,'(A)') ' added point '//TRIM(int2lstr(l))//' as fairlead for turbine '//trim(int2lstr(J)) + end if else @@ -1212,6 +1288,16 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er !also set number of attached lines to zero initially m%PointList(l)%nAttached = 0 + ! write body information to log file + if (p%writeLog > 1) then + write(p%UnLog, '(A)' ) " - Point"//trim(num2lstr(l))//":" + write(p%UnLog, '(A12,I2)' ) " id : ", m%PointList(l)%IdNum + write(p%UnLog, '(A12,I2)' ) " type: ", m%PointList(l)%typeNum + write(p%UnLog, '(A12,f12.4)') " v : ", m%PointList(l)%pointV + write(p%UnLog, '(A12,f12.4)') " m : ", m%PointList(l)%pointM + write(p%UnLog, '(A12,f12.4)') " CdA : ", m%PointList(l)%pointCdA + write(p%UnLog, '(A12,f12.4)') " Ca : ", m%PointList(l)%pointCa + end if ! check for sequential IdNums IF ( m%PointList(l)%IdNum .NE. l ) THEN @@ -1390,6 +1476,15 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! specify IdNum of line for error checking m%LineList(l)%IdNum = l + if (p%writeLog > 1) then + write(p%UnLog, '(A)' ) " - Line"//trim(num2lstr(m%LineList(l)%IdNum))//":" + write(p%UnLog, '(A15,I2)' ) " ID : ", m%LineList(l)%IdNum + write(p%UnLog, '(A15,A)' ) " Type : ", trim(m%LineTypeList(m%LineList(l)%PropsIdNum)%name) + write(p%UnLog, '(A15,f12.4)') " Len : ", m%LineList(l)%UnstrLen + write(p%UnLog, '(A15,A)' ) " Node A : ", " "//tempString2 + write(p%UnLog, '(A15,A)' ) " Node B : ", " "//tempString3 + write(p%UnLog, '(A15,I2)' ) " NumSegs: ", m%LineList(l)%N + end if ! check for sequential IdNums IF ( m%LineList(l)%IdNum .NE. l ) THEN @@ -1443,11 +1538,20 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er if (m%LineList( TempIDnums(J) )%CtrlChan == 0) then ! ensure line doesn't already have a CtrlChan assigned m%LineList( TempIDnums(J) )%CtrlChan = Itemp CALL WrScr('Assigned Line '//TRIM(Int2LStr(TempIDnums(J)))//' to control channel '//TRIM(Int2LStr(Itemp))) + if (p%writeLog > 0) then + write(p%UnLog,'(A)') 'Assigned Line '//TRIM(Int2LStr(TempIDnums(J)))//' to control channel '//TRIM(Int2LStr(Itemp)) + end if else CALL WrScr('Error: Line '//TRIM(Int2LStr(TempIDnums(J)))//' already is assigned to control channel '//TRIM(Int2LStr(m%LineList( TempIDnums(J) )%CtrlChan))//' so cannot also be assigned to channel '//TRIM(Int2LStr(Itemp))) + if (p%writeLog > 0) then + write(p%UnLog,'(A)') 'Error: Line '//TRIM(Int2LStr(TempIDnums(J)))//' already is assigned to control channel '//TRIM(Int2LStr(m%LineList( TempIDnums(J) )%CtrlChan))//' so cannot also be assigned to channel '//TRIM(Int2LStr(Itemp)) + end if end if else CALL WrScr('Error: Line ID '//TRIM(Int2LStr(TempIDnums(J)))//' of CtrlChan '//TRIM(Int2LStr(Itemp))//' is out of range') + if (p%writeLog > 0) then + write(p%UnLog,'(A)') 'Error: Line ID '//TRIM(Int2LStr(TempIDnums(J)))//' of CtrlChan '//TRIM(Int2LStr(Itemp))//' is out of range' + end if end if END DO @@ -1459,6 +1563,9 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er else if (INDEX(Line, "FAILURE") > 0) then ! if failure conditions header CALL WrScr(" Warning: Failure capabilities are not yet implemented in MoorDyn.") + if (p%writeLog > 0) then + write(p%UnLog,'(A)') " Warning: Failure capabilities are not yet implemented in MoorDyn." + end if ! skip following two lines (label line and unit line) Line = NextLine(i) @@ -1534,7 +1641,12 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er CALL CheckError( ErrStat2, ErrMsg2 ) IF (ErrStat >= AbortErrLev) RETURN - + if (p%writeLog > 1) then + write(p%UnLog, '(A)' ) " - Outputs List:" + DO J = 1, p%NumOuts + write(p%UnLog, '(A)' ) " "//OutList(J) + END DO + end if !------------------------------------------------------------------------------------------- else ! otherwise ignore this line that isn't a recognized header line and read the next line Line = NextLine(i) @@ -1573,7 +1685,10 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! p%NAnchs = 0 ! this is the number of "fixed" type Points. <<<<<<<<<<<<<< CALL WrScr(trim(Num2LStr(p%nLines))//' lines, '//trim(Num2LStr(p%NPoints))//' points, '//trim(Num2LStr(p%nRods))//' rods, '//trim(Num2LStr(p%nBodies))//' bodies.') - + if (p%writeLog > 0) then + write(p%UnLog, '(A)') NewLine + write(p%UnLog, '(A)') ' Created mooring system: '//trim(Num2LStr(p%nLines))//' lines, '//trim(Num2LStr(p%NPoints))//' points, '//trim(Num2LStr(p%nRods))//' rods, '//trim(Num2LStr(p%nBodies))//' bodies.' + end if @@ -1614,7 +1729,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! write system description to log file if (p%writeLog > 1) then - write(p%UnLog, '(A)') "----- MoorDyn Model Summary (to be written) -----" + write(p%UnLog, '(A)') "----- MoorDyn Model Summary (unfinished) -----" end if @@ -1763,20 +1878,20 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er J = J + 1 rRef = m%BodyList(m%CpldBodyIs(l,iTurb))%r6 ! for now set reference position as per input file <<< - !OrMatRef = CALL MeshPositionNode(u%CoupledKinematics(iTurb), J, rRef(1:3), ErrStat2, ErrMsg2) ! defaults to identity orientation matrix - !TODO: >>> should also maybe set reference orientation (which might make part of a couple lines down redundant) <<< - - ! calculate initial point relative position, adjusted due to initial platform translations - u%CoupledKinematics(iTurb)%TranslationDisp(:,J) = InitInp%PtfmInit(1:3,iTurb) - rRef(1:3) - - OrMat2 = MATMUL(OrMat, TRANSPOSE( EulerConstruct( rRef(4:6)))) ! combine the Body's relative orientation with the turbine's initial orientation - u%CoupledKinematics(iTurb)%Orientation(:,:,J) = OrMat2 ! set the result as the current orientation of the body <<< ! set absolute initial positions in MoorDyn - m%BodyList(m%CpldBodyIs(l,iTurb))%r6(1:3) = u%CoupledKinematics(iTurb)%Position(:,J) + u%CoupledKinematics(iTurb)%TranslationDisp(:,J) + p%TurbineRefPos(:,iTurb) - m%BodyList(m%CpldBodyIs(l,iTurb))%r6(4:6) = EulerExtract(MATMUL(OrMat, TRANSPOSE( EulerConstruct( rRef(4:6))))) ! apply rotation from PtfmInit onto input file's body orientation to get its true initial orientation + IF (p%Standalone /= 1) THEN + !TODO: >>> should also maybe set reference orientation (which might make part of a couple lines down redundant) <<< + OrMat2 = MATMUL(OrMat, ( EulerConstruct( rRef(4:6)))) ! combine the Body's relative orientation with the turbine's initial orientation + u%CoupledKinematics(iTurb)%Orientation(:,:,J) = OrMat2 ! set the result as the current orientation of the body <<< + + ! calculate initial point relative position, adjusted due to initial platform translations + u%CoupledKinematics(iTurb)%TranslationDisp(:,J) = InitInp%PtfmInit(1:3,iTurb) - rRef(1:3) + m%BodyList(m%CpldBodyIs(l,iTurb))%r6(1:3) = u%CoupledKinematics(iTurb)%Position(:,J) + u%CoupledKinematics(iTurb)%TranslationDisp(:,J) + p%TurbineRefPos(:,iTurb) + m%BodyList(m%CpldBodyIs(l,iTurb))%r6(4:6) = EulerExtract(OrMat2) ! apply rotation from PtfmInit onto input file's body orientation to get its true initial orientation + ENDIF CALL MeshConstructElement(u%CoupledKinematics(iTurb), ELEMENT_POINT, ErrStat2, ErrMsg2, J) ! set node as point element @@ -1789,23 +1904,24 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er J = J + 1 rRef = m%RodList(m%CpldRodIs(l,iTurb))%r6 ! for now set reference position as per input file <<< - OrMatRef = TRANSPOSE( m%RodList(m%CpldRodIs(l,iTurb))%OrMat ) ! for now set reference orientation as per input file <<< - CALL MeshPositionNode(u%CoupledKinematics(iTurb), J, rRef(1:3), ErrStat2, ErrMsg2, OrMatRef) ! assign the reference position and orientation - - ! calculate initial point relative position, adjusted due to initial platform rotations and translations <<< could convert to array math - u%CoupledKinematics(iTurb)%TranslationDisp(1,J) = InitInp%PtfmInit(1,iTurb) + OrMat(1,1)*rRef(1) + OrMat(2,1)*rRef(2) + OrMat(3,1)*rRef(3) - rRef(1) - u%CoupledKinematics(iTurb)%TranslationDisp(2,J) = InitInp%PtfmInit(2,iTurb) + OrMat(1,2)*rRef(1) + OrMat(2,2)*rRef(2) + OrMat(3,2)*rRef(3) - rRef(2) - u%CoupledKinematics(iTurb)%TranslationDisp(3,J) = InitInp%PtfmInit(3,iTurb) + OrMat(1,3)*rRef(1) + OrMat(2,3)*rRef(2) + OrMat(3,3)*rRef(3) - rRef(3) - - OrMat2 = MATMUL(OrMat, TRANSPOSE( EulerConstruct( rRef(4:6)))) ! combine the Rod's relative orientation with the turbine's initial orientation - u%CoupledKinematics(iTurb)%Orientation(:,:,J) = OrMat2 ! set the result as the current orientation of the rod <<< - + ! set absolute initial positions in MoorDyn - m%RodList(m%CpldRodIs(l,iTurb))%r6(1:3) = u%CoupledKinematics(iTurb)%Position(:,J) + u%CoupledKinematics(iTurb)%TranslationDisp(:,J) + p%TurbineRefPos(:,iTurb) - m%RodList(m%CpldRodIs(l,iTurb))%r6(4:6) = EulerExtract(MATMUL(OrMat, OrMatRef)) ! apply rotation from PtfmInit onto input file's rod orientation to get its true initial orientation + IF (p%Standalone /= 1) THEN + OrMatRef = ( m%RodList(m%CpldRodIs(l,iTurb))%OrMat ) ! for now set reference orientation as per input file <<< + CALL MeshPositionNode(u%CoupledKinematics(iTurb), J, rRef(1:3), ErrStat2, ErrMsg2, OrMatRef) ! assign the reference position and orientation + OrMat2 = MATMUL(OrMat, OrMatRef) ! combine the Rod's relative orientation with the turbine's initial orientation + u%CoupledKinematics(iTurb)%Orientation(:,:,J) = OrMat2 ! set the result as the current orientation of the rod <<< + + ! calculate initial point relative position, adjusted due to initial platform rotations and translations <<< could convert to array math + u%CoupledKinematics(iTurb)%TranslationDisp(1,J) = InitInp%PtfmInit(1,iTurb) + OrMat(1,1)*rRef(1) + OrMat(2,1)*rRef(2) + OrMat(3,1)*rRef(3) - rRef(1) + u%CoupledKinematics(iTurb)%TranslationDisp(2,J) = InitInp%PtfmInit(2,iTurb) + OrMat(1,2)*rRef(1) + OrMat(2,2)*rRef(2) + OrMat(3,2)*rRef(3) - rRef(2) + u%CoupledKinematics(iTurb)%TranslationDisp(3,J) = InitInp%PtfmInit(3,iTurb) + OrMat(1,3)*rRef(1) + OrMat(2,3)*rRef(2) + OrMat(3,3)*rRef(3) - rRef(3) + m%RodList(m%CpldRodIs(l,iTurb))%r6(1:3) = u%CoupledKinematics(iTurb)%Position(:,J) + u%CoupledKinematics(iTurb)%TranslationDisp(:,J) + p%TurbineRefPos(:,iTurb) + m%RodList(m%CpldRodIs(l,iTurb))%r6(4:6) = MATMUL(OrMat2 , (/0.0, 0.0, 1.0/) ) ! apply rotation from PtfmInit onto input file's rod orientation to get its true initial orientation + ENDIF ! >>> still need to set Rod initial orientations accounting for PtfmInit rotation <<< - + CALL MeshConstructElement(u%CoupledKinematics(iTurb), ELEMENT_POINT, ErrStat2, ErrMsg2, J) ! lastly, do this to set the attached line endpoint positions: @@ -1817,16 +1933,16 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! set reference position as per input file <<< what about turbine positions in array? rRef(1:3) = m%PointList(m%CpldPointIs(l,iTurb))%r - CALL MeshPositionNode(u%CoupledKinematics(iTurb), J, rRef(1:3), ErrStat2, ErrMsg2) - - ! calculate initial point relative position, adjusted due to initial platform rotations and translations <<< could convert to array math - u%CoupledKinematics(iTurb)%TranslationDisp(1,J) = InitInp%PtfmInit(1,iTurb) + OrMat(1,1)*rRef(1) + OrMat(2,1)*rRef(2) + OrMat(3,1)*rRef(3) - rRef(1) - u%CoupledKinematics(iTurb)%TranslationDisp(2,J) = InitInp%PtfmInit(2,iTurb) + OrMat(1,2)*rRef(1) + OrMat(2,2)*rRef(2) + OrMat(3,2)*rRef(3) - rRef(2) - u%CoupledKinematics(iTurb)%TranslationDisp(3,J) = InitInp%PtfmInit(3,iTurb) + OrMat(1,3)*rRef(1) + OrMat(2,3)*rRef(2) + OrMat(3,3)*rRef(3) - rRef(3) - + ! set absolute initial positions in MoorDyn - m%PointList(m%CpldPointIs(l,iTurb))%r = u%CoupledKinematics(iTurb)%Position(:,J) + u%CoupledKinematics(iTurb)%TranslationDisp(:,J) + p%TurbineRefPos(:,iTurb) - + IF (p%Standalone /= 1) THEN + CALL MeshPositionNode(u%CoupledKinematics(iTurb), J, rRef(1:3), ErrStat2, ErrMsg2) + ! calculate initial point relative position, adjusted due to initial platform rotations and translations <<< could convert to array math + u%CoupledKinematics(iTurb)%TranslationDisp(1,J) = InitInp%PtfmInit(1,iTurb) + OrMat(1,1)*rRef(1) + OrMat(2,1)*rRef(2) + OrMat(3,1)*rRef(3) - rRef(1) + u%CoupledKinematics(iTurb)%TranslationDisp(2,J) = InitInp%PtfmInit(2,iTurb) + OrMat(1,2)*rRef(1) + OrMat(2,2)*rRef(2) + OrMat(3,2)*rRef(3) - rRef(2) + u%CoupledKinematics(iTurb)%TranslationDisp(3,J) = InitInp%PtfmInit(3,iTurb) + OrMat(1,3)*rRef(1) + OrMat(2,3)*rRef(2) + OrMat(3,3)*rRef(3) - rRef(3) + m%PointList(m%CpldPointIs(l,iTurb))%r = u%CoupledKinematics(iTurb)%Position(:,J) + u%CoupledKinematics(iTurb)%TranslationDisp(:,J) + p%TurbineRefPos(:,iTurb) + ENDIF CALL MeshConstructElement(u%CoupledKinematics(iTurb), ELEMENT_POINT, ErrStat2, ErrMsg2, J) ! lastly, do this to set the attached line endpoint positions: @@ -1942,7 +2058,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! >>> maybe this should be skipped <<<< - ! Go through Bodys and write the coordinates to the state vector + ! Go through free Bodys (including pinned) and write the coordinates to the state vector DO l = 1,p%nFreeBodies CALL Body_Initialize(m%BodyList(m%FreeBodyIs(l)), x%states(m%BodyStateIs1(l) : m%BodyStateIsN(l)), m) END DO @@ -2035,7 +2151,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! if log file, compute and write some object properties ! ------------------------------------------------------------------- if (p%writeLog > 1) then - + write(p%UnLog, '(A)' ) "Values after intialization before dynamic relaxation" write(p%UnLog, '(A)' ) " Bodies:" DO l = 1,p%nBodies write(p%UnLog, '(A)' ) " Body"//trim(num2lstr(l))//":" @@ -2045,21 +2161,21 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er write(p%UnLog, '(A)' ) " Rods:" DO l = 1,p%nRods write(p%UnLog, '(A)' ) " Rod"//trim(num2lstr(l))//":" - ! m%RodList(l) + write(p%UnLog, '(A12, f12.4)') " mass: ", m%RodList(l)%M6net(1,1) + write(p%UnLog, '(A17, A)') " direction: ", trim(num2lstr(m%RodList(l)%q(1)))//", "//trim(num2lstr(m%RodList(l)%q(2)))//", "//trim(num2lstr(m%RodList(l)%q(3))) END DO write(p%UnLog, '(A)' ) " Points:" DO l = 1,p%nFreePoints write(p%UnLog, '(A)' ) " Point"//trim(num2lstr(l))//":" - ! m%PointList(l) + write(p%UnLog, '(A12, f12.4)') " mass: ", m%PointList(l)%M END DO write(p%UnLog, '(A)' ) " Lines:" DO l = 1,p%nLines write(p%UnLog, '(A)' ) " Line"//trim(num2lstr(l))//":" - ! m%LineList(l) END DO - + write(p%UnLog, '(A)') "--------- End of Model Summary --------- "//NewLine end if @@ -2071,6 +2187,9 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er if (InputFileDat%TMaxIC > 0.0_DbKi) then CALL WrScr(" Finalizing initial conditions using dynamic relaxation."//NewLine) ! newline because next line writes over itself + if (p%writeLog > 0) then + write(p%UnLog,'(A)') "Finalizing initial conditions using dynamic relaxation."//NewLine + end if ! boost drag coefficient of each line type <<<<<<<< does this actually do anything or do lines hold these coefficients??? DO I = 1, p%nLines @@ -2139,6 +2258,10 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er IF (ErrStat == ErrID_Fatal) THEN CALL WrScr("NaN detected at time "//TRIM(Num2LStr(t))//" during MoorDyn's dynamic relaxation process.") + if (p%writeLog > 0) then + write(p%UnLog,'(A)') "NaN detected at time "//TRIM(Num2LStr(t))//" during MoorDyn's dynamic relaxation process."//NewLine + end if + IF (wordy > 1) THEN print *, "Here is the state vector: " print *, x%states @@ -2190,11 +2313,18 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er END DO IF (Converged == 1) THEN ! if we made it with all cases satisfying the threshold - CALL WrScr('') ! serves as line break from write over command in previous printed line CALL WrScr(' Fairlead tensions converged to '//trim(Num2LStr(100.0*InputFileDat%threshIC))//'% after '//trim(Num2LStr(t))//' seconds.') + if (p%writeLog > 0) then + write(p%UnLog,'(A)') '' + write(p%UnLog,'(A)') ' Fairlead tensions converged to '//trim(Num2LStr(100.0*InputFileDat%threshIC))//'% after '//trim(Num2LStr(t))//' seconds.'//NewLine + end if DO l = 1, p%nLines CALL WrScr(' Fairlead tension: '//trim(Num2LStr(FairTensIC(l,1)))) CALL WrScr(' Fairlead forces: '//trim(Num2LStr(m%LineList(l)%Fnet(1, m%LineList(l)%N)))//', '//trim(Num2LStr(m%LineList(l)%Fnet(2, m%LineList(l)%N)))//', '//trim(Num2LStr(m%LineList(l)%Fnet(3, m%LineList(l)%N)))) + if (p%writeLog > 0) then + write(p%UnLog,'(A)') ' Fairlead tension: '//trim(Num2LStr(FairTensIC(l,1))) + write(p%UnLog,'(A)') ' Fairlead forces: '//trim(Num2LStr(m%LineList(l)%Fnet(1, m%LineList(l)%N)))//', '//trim(Num2LStr(m%LineList(l)%Fnet(2, m%LineList(l)%N)))//', '//trim(Num2LStr(m%LineList(l)%Fnet(3, m%LineList(l)%N))) + end if ENDDO EXIT ! break out of the time stepping loop END IF @@ -2203,6 +2333,11 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er IF (I == ceiling(InputFileDat%TMaxIC/InputFileDat%dtIC) ) THEN CALL WrScr('') ! serves as line break from write over command in previous printed line CALL WrScr(' Fairlead tensions did not converge within TMaxIC='//trim(Num2LStr(InputFileDat%TMaxIC))//' seconds.') + if (p%writeLog > 0) then + write(p%UnLog,'(A)') '' + write(p%UnLog,'(A)') ' Fairlead tensions did not converge within TMaxIC='//trim(Num2LStr(InputFileDat%TMaxIC))//' seconds.' + end if + !ErrStat = ErrID_Warn !ErrMsg = ' MD_Init: ran dynamic convergence to TMaxIC without convergence' END IF @@ -2247,6 +2382,14 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er endif CALL WrScr(' MoorDyn initialization completed.') + if (p%writeLog > 0) then + write(p%UnLog, '(A)') NewLine//"MoorDyn initialization completed."//NewLine + if (ErrStat /= ErrID_None) then + write(p%UnLog, '(A34)') "Initalization Errors and Warnings:" + write(p%UnLog, '(A)' ) ErrMsg + end if + write(p%UnLog, '(A)') NewLine + end if m%LastOutTime = -1.0_DbKi ! set to nonzero to ensure that output happens at the start of simulation at t=0 @@ -2328,7 +2471,6 @@ END SUBROUTINE CheckError SUBROUTINE CleanUp() ! ErrStat = ErrID_Fatal call MD_DestroyInputFileType( InputFileDat, ErrStat2, ErrMsg2 ) ! Ignore any error messages from this - IF (p%UnLog > 0_IntKi) CLOSE( p%UnLog ) ! Remove this when the log file is kept open during the full simulation END SUBROUTINE !> If for some reason the file is truncated, it is possible to get into an infinite loop @@ -2440,6 +2582,9 @@ SUBROUTINE MD_UpdateStates( t, n, u, t_array, p, x, xd, z, other, m, ErrStat, Er IF (ErrStat == ErrID_Fatal) THEN CALL WrScr("NaN detected at time "//TRIM(Num2LStr(t2))//" in MoorDyn.") + if (p%writeLog > 0) then + write(p%UnLog,'(A)') "NaN detected at time "//TRIM(Num2LStr(t2))//" in MoorDyn." + end if IF (wordy > 1) THEN print *, ". Here is the state vector: " print *, x%states @@ -2471,6 +2616,9 @@ SUBROUTINE MD_UpdateStates( t, n, u, t_array, p, x, xd, z, other, m, ErrStat, Er IF (ErrStat == ErrID_Fatal) THEN CALL WrScr("NaN detected at time "//TRIM(Num2LStr(t2))//" in MoorDyn.") + if (p%writeLog > 0) then + write(p%UnLog,'(A)') "NaN detected at time "//TRIM(Num2LStr(t2))//" in MoorDyn." + end if IF (wordy > 1) THEN print *, ". Here is the state vector: " print *, x%states @@ -2493,6 +2641,9 @@ SUBROUTINE CheckError(ErrId, Msg) ErrStat = MAX(ErrStat, ErrID) CALL WrScr( ErrMsg ) ! do this always or only if warning level? + if (p%writeLog > 0) then + write(p%UnLog,'(A)') ErrMsg + end if IF( ErrStat > ErrID_Warn ) THEN ! CALL MD_DestroyInput( u_interp, ErrStat, ErrMsg ) @@ -2697,6 +2848,9 @@ SUBROUTINE CheckError(ErrId, Msg) ErrStat = MAX(ErrStat, ErrID) CALL WrScr( ErrMsg ) ! do this always or only if warning level? <<<<<<<<<<<<<<<<<<<<<< probably should remove all instances + if (p%writeLog > 0) then + write(p%UnLog,'(A)') ErrMsg + end if ! IF( ErrStat > ErrID_Warn ) THEN ! CALL MD_DestroyContState( dxdt, ErrStat2, ErrMsg2) @@ -2849,6 +3003,9 @@ SUBROUTINE MD_CalcContStateDeriv( t, u, p, x, xd, z, other, m, dxdt, ErrStat, Er ErrStat = ErrID_Fatal ErrMsg = ' Active tension command will make a segment longer than the limit of twice its original length.' call WrScr(trim(Num2LStr(u%DeltaL(m%LineList(L)%CtrlChan)))//" is an increase of more than "//trim(Num2LStr(m%LineList(L)%UnstrLen / m%LineList(L)%N))) + if (p%writeLog > 0) then + write(p%UnLog,'(A)') trim(Num2LStr(u%DeltaL(m%LineList(L)%CtrlChan)))//" is an increase of more than "//trim(Num2LStr(m%LineList(L)%UnstrLen / m%LineList(L)%N)) + end if IF (wordy > 0) print *, u%DeltaL IF (wordy > 0) print*, m%LineList(L)%CtrlChan RETURN @@ -2857,6 +3014,9 @@ SUBROUTINE MD_CalcContStateDeriv( t, u, p, x, xd, z, other, m, dxdt, ErrStat, Er ErrStat = ErrID_Fatal ErrMsg = ' Active tension command will make a segment shorter than the limit of half its original length.' call WrScr(trim(Num2LStr(u%DeltaL(m%LineList(L)%CtrlChan)))//" is a reduction of more than half of "//trim(Num2LStr(m%LineList(L)%UnstrLen / m%LineList(L)%N))) + if (p%writeLog > 0) then + write(p%UnLog,'(A)') trim(Num2LStr(u%DeltaL(m%LineList(L)%CtrlChan)))//" is a reduction of more than half of "//trim(Num2LStr(m%LineList(L)%UnstrLen / m%LineList(L)%N)) + end if IF (wordy > 0) print *, u%DeltaL IF (wordy > 0) print*, m%LineList(L)%CtrlChan RETURN @@ -2965,19 +3125,23 @@ SUBROUTINE MD_CalcContStateDeriv( t, u, p, x, xd, z, other, m, dxdt, ErrStat, Er DO iTurb = 1,p%nTurbines DO l = 1,p%nCpldPoints(iTurb) - ! >>>>>>>> here we should pass along accelerations and include inertial loads in the calculation!!! <<>>>>>>> here we should pass along accelerations and include inertial loads in the calculation!!! << 0) then + write(p%UnLog,'(A)') ErrMsg + end if END IF @@ -3801,27 +3968,39 @@ SUBROUTINE Init_Jacobian_x() idx = 0 ! Free bodies DO l = 1,p%nFreeBodies ! Body m%BodyList(m%FreeBodyIs(l)) - p%dx(idx+1:idx+3) = dl_slack_min ! body displacement [m] - p%dx(idx+4:idx+6) = 0.02 ! body rotation [rad] - ! corresponds to state indices: (m%BodyStateIs1(l)+6:m%BodyStateIs1(l)+11) - InitOut%LinNames_x(idx+1) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' Px, m' - InitOut%LinNames_x(idx+2) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' Py, m' - InitOut%LinNames_x(idx+3) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' Pz, m' - InitOut%LinNames_x(idx+4) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' rot_x, rad' - InitOut%LinNames_x(idx+5) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' rot_y, rad' - InitOut%LinNames_x(idx+6) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' rot_z, rad' - p%dxIdx_map2_xStateIdx(idx+1) = m%BodyStateIs1(l)+6 ! x%state index for Px - p%dxIdx_map2_xStateIdx(idx+2) = m%BodyStateIs1(l)+7 ! x%state index for Py - p%dxIdx_map2_xStateIdx(idx+3) = m%BodyStateIs1(l)+8 ! x%state index for Pz - p%dxIdx_map2_xStateIdx(idx+4) = m%BodyStateIs1(l)+9 ! x%state index for rot_x - p%dxIdx_map2_xStateIdx(idx+5) = m%BodyStateIs1(l)+10 ! x%state index for rot_y - p%dxIdx_map2_xStateIdx(idx+6) = m%BodyStateIs1(l)+11 ! x%state index for rot_z - idx = idx + 6 + if (m%BodyList(m%FreeBodyIs(l))%typeNum == 2) then ! Coupled pinned body + p%dx(idx+4:idx+6) = 0.02 ! body rotation [rad] + ! corresponds to state indices: (m%BodyStateIs1(l)+6:m%BodyStateIs1(l)+8) + InitOut%LinNames_x(idx+1) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' rot_x, rad' + InitOut%LinNames_x(idx+2) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' rot_y, rad' + InitOut%LinNames_x(idx+3) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' rot_z, rad' + p%dxIdx_map2_xStateIdx(idx+4) = m%BodyStateIs1(l)+3 ! x%state index for rot_x + p%dxIdx_map2_xStateIdx(idx+5) = m%BodyStateIs1(l)+4 ! x%state index for rot_y + p%dxIdx_map2_xStateIdx(idx+6) = m%BodyStateIs1(l)+5 ! x%state index for rot_z + idx = idx + 3 + else ! free body + p%dx(idx+1:idx+3) = dl_slack_min ! body displacement [m] + p%dx(idx+4:idx+6) = 0.02 ! body rotation [rad] + ! corresponds to state indices: (m%BodyStateIs1(l)+6:m%BodyStateIs1(l)+11) + InitOut%LinNames_x(idx+1) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' Px, m' + InitOut%LinNames_x(idx+2) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' Py, m' + InitOut%LinNames_x(idx+3) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' Pz, m' + InitOut%LinNames_x(idx+4) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' rot_x, rad' + InitOut%LinNames_x(idx+5) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' rot_y, rad' + InitOut%LinNames_x(idx+6) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' rot_z, rad' + p%dxIdx_map2_xStateIdx(idx+1) = m%BodyStateIs1(l)+6 ! x%state index for Px + p%dxIdx_map2_xStateIdx(idx+2) = m%BodyStateIs1(l)+7 ! x%state index for Py + p%dxIdx_map2_xStateIdx(idx+3) = m%BodyStateIs1(l)+8 ! x%state index for Pz + p%dxIdx_map2_xStateIdx(idx+4) = m%BodyStateIs1(l)+9 ! x%state index for rot_x + p%dxIdx_map2_xStateIdx(idx+5) = m%BodyStateIs1(l)+10 ! x%state index for rot_y + p%dxIdx_map2_xStateIdx(idx+6) = m%BodyStateIs1(l)+11 ! x%state index for rot_z + idx = idx + 6 + endif END DO ! Rods DO l = 1,p%nFreeRods ! Rod m%RodList(m%FreeRodIs(l)) - if (m%RodList(m%FreeRodIs(l))%typeNum == 1) then ! pinned rod + if (abs(m%RodList(m%FreeRodIs(l))%typeNum) == 1) then ! pinned rod p%dx(idx+1:idx+3) = 0.02 ! rod rotation [rad] ! corresponds to state indices: (m%RodStateIs1(l)+3:m%RodStateIs1(l)+5) InitOut%LinNames_x(idx+1) = 'Rod '//trim(num2lstr(m%FreeRodIs(l)))//' rot_x, rad' @@ -3885,27 +4064,39 @@ SUBROUTINE Init_Jacobian_x() !----------------- ! Free bodies DO l = 1,p%nFreeBodies ! Body m%BodyList(m%FreeBodyIs(l)) - ! corresponds to state indices: (m%BodyStateIs1(l):m%BodyStateIs1(l)+5) - p%dx(idx+1:idx+3) = 0.1 ! body translational velocity [m/s] - p%dx(idx+4:idx+6) = 0.1 ! body rotational velocity [rad/s] - InitOut%LinNames_x(idx+1) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' Vx, m/s' - InitOut%LinNames_x(idx+2) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' Vy, m/s' - InitOut%LinNames_x(idx+3) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' Vz, m/s' - InitOut%LinNames_x(idx+4) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' omega_x, rad/s' - InitOut%LinNames_x(idx+5) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' omega_y, rad/s' - InitOut%LinNames_x(idx+6) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' omega_z, rad/s' - p%dxIdx_map2_xStateIdx(idx+1) = m%BodyStateIs1(l)+0 ! x%state index for Rx - p%dxIdx_map2_xStateIdx(idx+2) = m%BodyStateIs1(l)+1 ! x%state index for Ry - p%dxIdx_map2_xStateIdx(idx+3) = m%BodyStateIs1(l)+2 ! x%state index for Rz - p%dxIdx_map2_xStateIdx(idx+4) = m%BodyStateIs1(l)+3 ! x%state index for omega_x - p%dxIdx_map2_xStateIdx(idx+5) = m%BodyStateIs1(l)+4 ! x%state index for omega_y - p%dxIdx_map2_xStateIdx(idx+6) = m%BodyStateIs1(l)+5 ! x%state index for omega_z - idx = idx + 6 + if (m%BodyList(m%FreeBodyIs(l))%typeNum == 2) then ! Coupled pinned body + ! corresponds to state indices: (m%BodyStateIs1(l):m%BodyStateIs1(l)+5) + p%dx(idx+1:idx+3) = 0.1 ! body rotational velocity [rad/s] + InitOut%LinNames_x(idx+1) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' omega_x, rad/s' + InitOut%LinNames_x(idx+2) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' omega_y, rad/s' + InitOut%LinNames_x(idx+3) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' omega_z, rad/s' + p%dxIdx_map2_xStateIdx(idx+1) = m%BodyStateIs1(l)+0 ! x%state index for omega_x + p%dxIdx_map2_xStateIdx(idx+2) = m%BodyStateIs1(l)+1 ! x%state index for omega_y + p%dxIdx_map2_xStateIdx(idx+3) = m%BodyStateIs1(l)+2 ! x%state index for omega_z + idx = idx + 3 + else !Free body + ! corresponds to state indices: (m%BodyStateIs1(l):m%BodyStateIs1(l)+5) + p%dx(idx+1:idx+3) = 0.1 ! body translational velocity [m/s] + p%dx(idx+4:idx+6) = 0.1 ! body rotational velocity [rad/s] + InitOut%LinNames_x(idx+1) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' Vx, m/s' + InitOut%LinNames_x(idx+2) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' Vy, m/s' + InitOut%LinNames_x(idx+3) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' Vz, m/s' + InitOut%LinNames_x(idx+4) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' omega_x, rad/s' + InitOut%LinNames_x(idx+5) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' omega_y, rad/s' + InitOut%LinNames_x(idx+6) = 'Body '//trim(num2lstr(m%FreeBodyIs(l)))//' omega_z, rad/s' + p%dxIdx_map2_xStateIdx(idx+1) = m%BodyStateIs1(l)+0 ! x%state index for Rx + p%dxIdx_map2_xStateIdx(idx+2) = m%BodyStateIs1(l)+1 ! x%state index for Ry + p%dxIdx_map2_xStateIdx(idx+3) = m%BodyStateIs1(l)+2 ! x%state index for Rz + p%dxIdx_map2_xStateIdx(idx+4) = m%BodyStateIs1(l)+3 ! x%state index for omega_x + p%dxIdx_map2_xStateIdx(idx+5) = m%BodyStateIs1(l)+4 ! x%state index for omega_y + p%dxIdx_map2_xStateIdx(idx+6) = m%BodyStateIs1(l)+5 ! x%state index for omega_z + idx = idx + 6 + endif END DO ! Rods DO l = 1,p%nFreeRods ! Rod m%RodList(m%FreeRodIs(l)) - if (m%RodList(m%FreeRodIs(l))%typeNum == 1) then ! pinned rod + if (abs(m%RodList(m%FreeRodIs(l))%typeNum) == 1) then ! pinned rod ! corresponds to state indices: (m%RodStateIs1(l):m%RodStateIs1(l)+2) p%dx(idx+1:idx+3) = 0.1 ! body rotational velocity [rad/s] InitOut%LinNames_x(idx+1) = 'Rod '//trim(num2lstr(m%FreeRodIs(l)))//' omega_x, rad/s' diff --git a/modules/moordyn/src/MoorDyn_Body.f90 b/modules/moordyn/src/MoorDyn_Body.f90 index 74fe1ec60a..acf6f92098 100644 --- a/modules/moordyn/src/MoorDyn_Body.f90 +++ b/modules/moordyn/src/MoorDyn_Body.f90 @@ -155,12 +155,25 @@ SUBROUTINE Body_Initialize(Body, states, m) INTEGER(IntKi) :: l ! index of segments or nodes along line REAL(DbKi) :: dummyStates(12) ! dummy vector to mimic states when initializing a rigidly attached rod - - ! assign initial body kinematics to state vector - states(7:12) = Body%r6 - states(1:6 ) = Body%v6 - + IF (wordy > 0) print *, "initializing Body ", Body%idNum + + ! the r6 and v6 vectors should have already been set + ! r and rd of ends have already been set by setup function or by parent object <<<<< right? <<<<< + + if (Body%typeNum == 0) then ! free body type + + ! assign initial body kinematics to state vector + states(1:6 ) = Body%v6 ! zero velocities for initialization (set to 0 in Body_Setup) + states(7:12) = Body%r6 + + else if (Body%typeNum ==2 ) then ! pinned rod type (coupled or attached to something previously via setPinKin) + + states(1:3) = Body%v6(4:6) ! zero velocities for initialization (set to 0 in Body_Setup) + states(4:6) = Body%r6(4:6) ! body orentations + + end if + ! set positions of any dependent points and rods now (before they are initialized) CALL Body_SetDependentKin(Body, 0.0_DbKi, m) @@ -203,12 +216,12 @@ END SUBROUTINE Body_InitializeUnfree ! set kinematics for Bodies if they are coupled (or ground) !-------------------------------------------------------------- - SUBROUTINE Body_SetKinematics(Body, r_in, v_in, a_in, t, m) + SUBROUTINE Body_SetKinematics(Body, r6_in, v6_in, a6_in, t, m) Type(MD_Body), INTENT(INOUT) :: Body ! the Body object - Real(DbKi), INTENT(IN ) :: r_in(6) ! 6-DOF position - Real(DbKi), INTENT(IN ) :: v_in(6) ! 6-DOF velocity - Real(DbKi), INTENT(IN ) :: a_in(6) ! 6-DOF acceleration (only used for coupled rods) + Real(DbKi), INTENT(IN ) :: r6_in(6) ! 6-DOF position + Real(DbKi), INTENT(IN ) :: v6_in(6) ! 6-DOF velocity + Real(DbKi), INTENT(IN ) :: a6_in(6) ! 6-DOF acceleration Real(DbKi), INTENT(IN ) :: t ! instantaneous time TYPE(MD_MiscVarType), INTENT(INOUT) :: m ! passing along all mooring objects (for simplicity, since Bodies deal with Rods and Points) @@ -218,26 +231,24 @@ SUBROUTINE Body_SetKinematics(Body, r_in, v_in, a_in, t, m) ! store current time Body%time = t - ! if (abs(Body%typeNum) == 2) then ! body coupled in 6 DOF, or ground - Body%r6 = r_in - Body%v6 = v_in - Body%a6 = a_in + if (Body%typeNum == 2) then ! body pinned to coupling point + + ! set Body translational kinematics based on BCs (linear model for now) + Body%r6(1:3) = r6_in(1:3) + Body%v6(1:3) = v6_in(1:3) + Body%a6(1:3) = a6_in(1:3) + + ! Body rotations are left alone and will be handled, along with passing kinematics to dependent objects, by separate call to setState + + else ! body rigidly coupled to coupling point + Body%r6 = r6_in + Body%v6 = v6_in + Body%a6 = a6_in ! since this body has no states and all DOFs have been set, pass its kinematics to dependent attachments CALL Body_SetDependentKin(Body, t, m) - - ! else if (abs(Body%typeNum) == 1) then ! body pinned at reference point - ! - ! ! set Body *end A only* kinematics based on BCs (linear model for now) - ! Body%r6(1:3) = r_in(1:3) - ! Body%v6(1:3) = v_in(1:3) - ! - ! ! Body is pinned so only ref point posiiton is specified, rotations are left alone and will be - ! ! handled, along with passing kinematics to attached objects, by separate call to setState - ! - ! else - ! print *, "Error: Body_SetKinematics called for a free Body." ! <<< - ! end if + + end if END SUBROUTINE Body_SetKinematics !-------------------------------------------------------------- @@ -257,14 +268,26 @@ SUBROUTINE Body_SetState(Body, X, t, m) ! store current time Body%time = t + if (Body%typeNum == 0) then ! free Body type - - Body%r6 = X(7:12) ! get positions - Body%v6 = X(1:6) ! get velocities - + Body%r6 = X(7:12) ! get positions + Body%v6 = X(1:6) ! get velocities + + ! set positions of any dependent points and rods + CALL Body_SetDependentKin(Body, t, m) + + else if (Body%typeNum == 2) then + + Body%r6(4:6) = X(4:6) ! get positions + Body%v6(4:6) = X(1:3) ! get velocities - ! set positions of any dependent points and rods - CALL Body_SetDependentKin(Body, t, m) + + ! set positions of any dependent points and rods + CALL Body_SetDependentKin(Body, t, m) + + else + print *, "Error: Body::setState called for a non-free Body type in MoorDyn" ! <<< + end if END SUBROUTINE Body_SetState !-------------------------------------------------------------- @@ -336,6 +359,8 @@ SUBROUTINE Body_GetStateDeriv(Body, Xd, m, p) INTEGER(IntKi) :: J ! index + Real(DbKi) :: Fnet (6) ! net force and moment about reference point + Real(DbKi) :: acc(6) ! 6DOF acceleration vector Real(DbKi) :: y_temp (6) ! temporary vector for LU decomposition @@ -349,15 +374,35 @@ SUBROUTINE Body_GetStateDeriv(Body, Xd, m, p) CALL Body_DoRHS(Body, m, p) - ! solve for accelerations in [M]{a}={f} using LU decomposition - CALL LUsolve(6, Body%M, LU_temp, Body%F6net, y_temp, acc) + IF (Body%typeNum == 0) THEN ! Free body - ! fill in state derivatives - Xd(7:12) = Body%v6 ! dxdt = V (velocities) - Xd(1:6) = acc ! dVdt = a (accelerations) + ! solve for accelerations in [M]{a}={f} using LU decomposition + CALL LUsolve(6, Body%M, LU_temp, Body%F6net, y_temp, acc) - ! store accelerations in case they're useful as output - Body%a6 = acc + ! fill in state derivatives + Xd(7:12) = Body%v6 ! dxdt = V (velocities) + Xd(1:6) = acc ! dVdt = a (accelerations) + + ! store accelerations in case they're useful as output + Body%a6 = acc + + ELSE ! Pinned Body, 3 states (rotational only) + + ! Account for moment response due to inertial coupling + Fnet = Body%F6net + Fnet(4:6) = Fnet(4:6) - MATMUL(Body%M(4:6,1:3), Body%a6(1:3)) + + ! solve for accelerations in [M]{a}={f} using LU decomposition + CALL LUsolve(3, Body%M(4:6,4:6), LU_temp(4:6,4:6), Fnet(4:6), y_temp(4:6), acc(4:6)) + + ! fill in state derivatives + Xd(4:6) = Body%v6(4:6) ! dxdt = V (velocities) + Xd(1:3) = acc(4:6) ! dVdt = a (accelerations) + + ! store accelerations in case they're useful as output + Body%a6(4:6) = acc(4:6) + + ENDIF ! check for NaNs (should check all state derivatives, not just first 6) DO J = 1, 6 @@ -477,9 +522,29 @@ SUBROUTINE Body_GetCoupledForce(Body, Fnet_out, m, p) ! add inertial loads as appropriate if (Body%typeNum == -1) then - F6_iner = 0.0_DbKi !-MATMUL(Body%M, Body%a6) <<<<<<<< why does including F6_iner cause instability??? - Fnet_out = Body%F6net + F6_iner ! add inertial loads + if (p%inertialF == 1) then ! include inertial components + F6_iner = -MATMUL(Body%M, Body%a6) ! unstable in OpenFAST v4 and below becasue of loose coupling with ED and SD. Transients in acceleration can cause issues + else + ! When OpenFAST v5 is released w/ tight coupling, remove this hack and just use the inertial term above + F6_iner = 0.0 + endif + + Body%F6net = Body%F6net + F6_iner ! add inertial loads + Fnet_out = Body%F6net + + else if (Body%typeNum == 2) then ! pinned coupled body + + if (p%inertialF == 1) then ! include inertial components + ! inertial loads ... from input translational ... and solved rotational ... acceleration + F6_iner(1:3) = -MATMUL(Body%M(1:3,1:3), Body%a6(1:3)) - MATMUL(Body%M(1:3,4:6), Body%a6(4:6)) + else + F6_iner(1:3) = 0.0 + endif + Body%F6net(1:3) = Body%F6net(1:3) + F6_iner(1:3) ! add translational inertial loads + Body%F6net(4:6) = 0.0_DbKi + Fnet_out = Body%F6net + else print *, "ERROR, Body_GetCoupledForce called for wrong (non-coupled) body type in MoorDyn!" end if @@ -505,7 +570,7 @@ SUBROUTINE Body_AddPoint(Body, pointID, coords) Body%AttachedC(Body%nAttachedC) = pointID Body%rPointRel(:,Body%nAttachedC) = coords ! store relative position of point on body ELSE - Print*, "too many Points attached to Body ", Body%IdNum, " in MoorDyn!" + call WrScr("too many Points attached to Body "//trim(num2lstr(Body%IdNum))//" in MoorDyn!") END IF END SUBROUTINE Body_AddPoint @@ -536,7 +601,7 @@ SUBROUTINE Body_AddRod(Body, rodID, coords) Body%r6RodRel(4:6, Body%nAttachedR) = tempUnitVec ELSE - Print*, "too many rods attached to Body ", Body%IdNum, " in MoorDyn" + call WrScr("too many rods attached to Body "//trim(num2lstr(Body%IdNum))//" in MoorDyn") END IF END SUBROUTINE Body_AddRod diff --git a/modules/moordyn/src/MoorDyn_Driver.f90 b/modules/moordyn/src/MoorDyn_Driver.f90 index 2f76811f9e..3fbe8b580a 100644 --- a/modules/moordyn/src/MoorDyn_Driver.f90 +++ b/modules/moordyn/src/MoorDyn_Driver.f90 @@ -28,7 +28,7 @@ PROGRAM MoorDyn_Driver IMPLICIT NONE TYPE MD_Drvr_InitInput - LOGICAL :: Echo + ! LOGICAL :: Echo REAL(DbKi) :: Gravity REAL(DbKi) :: rhoW REAL(DbKi) :: WtrDepth @@ -90,7 +90,7 @@ PROGRAM MoorDyn_Driver INTEGER(IntKi) :: nt ! number of coupling time steps to use in simulation REAL(DbKi) :: t ! current time (s) - REAL(DbKi) :: tMax ! sim end time (s) + REAL(DbKi) :: TMax ! sim end time (s) REAL(DbKi) :: dtC ! fixed/constant global time step REAL(DbKi) :: frac ! fraction used in interpolation @@ -120,7 +120,7 @@ PROGRAM MoorDyn_Driver ErrMsg = "" ErrStat = ErrID_None - UnEcho=-1 + UnEcho=-1 ! set to -1 as echo is no longer used by MD UnIn =-1 ! TODO: Sort out error handling (two sets of flags currently used) @@ -162,7 +162,7 @@ PROGRAM MoorDyn_Driver MD_InitInp%RootName = drvrInitInp%OutRootName MD_InitInp%UsePrimaryInputFile = .TRUE. !MD_InitInp%PassedPrimaryInputData = - MD_InitInp%Echo = drvrInitInp%Echo + ! MD_InitInp%Echo = drvrInitInp%Echo !MD_InitInp%OutList = <<<< never used? MD_InitInp%Linearize = .FALSE. @@ -175,7 +175,7 @@ PROGRAM MoorDyn_Driver if (drvrInitInp%FarmSize > 0) then ! Check if this MoorDyn instance is being run from FAST.Farm (indicated by FarmSize > 0) nTurbines = drvrInitInp%FarmSize - else ! FarmSize==0 indicates normal, FAST module mode + else ! FarmSize==0 indicates normal, FAST module mode; FarmSize<0 indicates standalone mode nTurbines = 1 ! if a regular FAST module mode, we treat it like a nTurbine=1 farm case end if @@ -300,7 +300,7 @@ PROGRAM MoorDyn_Driver ! specify stepping details - nt = tMax/dtC - 1 ! number of coupling time steps + nt = TMax/dtC - 1 ! number of coupling time steps ! allocate space for processed motion array @@ -447,11 +447,11 @@ PROGRAM MoorDyn_Driver else - nt = tMax/dtC - 1 ! number of coupling time steps + nt = TMax/dtC - 1 ! number of coupling time steps end if CALL WrScr(" ") - call WrScr("Tmax - "//trim(Num2LStr(tMax))//" and nt="//trim(Num2LStr(nt))) + call WrScr("Tmax - "//trim(Num2LStr(TMax))//" and nt="//trim(Num2LStr(nt))) CALL WrScr(" ") @@ -490,9 +490,14 @@ PROGRAM MoorDyn_Driver i = 1 ! read first timestep data K = 1 ! the index of the coupling points in the input mesh CoupledKinematics J = 1 ! the starting index of the relevant DOFs in the input array + + IF (MD_InitInp%FarmSize < 0) THEN + MD_p%TurbineRefPos(:,iTurb) = 0.0 + ENDIF + ! any coupled bodies (type -1) DO l = 1,MD_p%nCpldBodies(iTurb) - MD_u(1)%CoupledKinematics(iTurb)%TranslationDisp(:,K) = r_in(i, J:J+2) - MD_u(1)%CoupledKinematics(iTurb)%Position(:,K) - MD_p%TurbineRefPos(:,iTurb) + MD_u(1)%CoupledKinematics(iTurb)%TranslationDisp(:,K) = r_in(i, J:J+2) - MD_u(1)%CoupledKinematics(iTurb)%Position(:,K) - MD_p%TurbineRefPos(:,iTurb) MD_u(1)%CoupledKinematics(iTurb)%Orientation( :,:,K) = EulerConstruct( r_in(i, J+3:J+5) ) ! full Euler angle approach MD_u(1)%CoupledKinematics(iTurb)%TranslationVel( :,K) = rd_in(i, J:J+2) MD_u(1)%CoupledKinematics(iTurb)%RotationVel( :,K) = rd_in(i, J+3:J+5) @@ -505,8 +510,7 @@ PROGRAM MoorDyn_Driver ! any coupled rods (type -1 or -2) >>> need to make rotations ignored if it's a pinned rod <<< DO l = 1,MD_p%nCpldRods(iTurb) - - MD_u(1)%CoupledKinematics(iTurb)%TranslationDisp(:,K) = r_in(i, J:J+2) - MD_u(1)%CoupledKinematics(iTurb)%Position(:,K) - MD_p%TurbineRefPos(:,iTurb) + MD_u(1)%CoupledKinematics(iTurb)%TranslationDisp(:,K) = r_in(i, J:J+2) - MD_u(1)%CoupledKinematics(iTurb)%Position(:,K) - MD_p%TurbineRefPos(:,iTurb) MD_u(1)%CoupledKinematics(iTurb)%Orientation( :,:,K) = EulerConstruct( r_in(i, J+3:J+5) ) MD_u(1)%CoupledKinematics(iTurb)%TranslationVel( :,K) = rd_in(i, J:J+2) MD_u(1)%CoupledKinematics(iTurb)%RotationVel( :,K) = rd_in(i, J+3:J+5) @@ -519,8 +523,7 @@ PROGRAM MoorDyn_Driver ! any coupled points (type -1) DO l = 1, MD_p%nCpldPoints(iTurb) - - MD_u(1)%CoupledKinematics(iTurb)%TranslationDisp(:,K) = r_in(i, J:J+2) - MD_u(1)%CoupledKinematics(iTurb)%Position(:,K) - MD_p%TurbineRefPos(:,iTurb) + MD_u(1)%CoupledKinematics(iTurb)%TranslationDisp(:,K) = r_in(i, J:J+2) - MD_u(1)%CoupledKinematics(iTurb)%Position(:,K) - MD_p%TurbineRefPos(:,iTurb) MD_u(1)%CoupledKinematics(iTurb)%TranslationVel( :,K) = rd_in(i, J:J+2) MD_u(1)%CoupledKinematics(iTurb)%TranslationAcc( :,K) = 0.0_DbKi !rdd_in(i, J:J+2) @@ -553,7 +556,7 @@ PROGRAM MoorDyn_Driver call WrScr("Doing time marching now...") - CALL SimStatus_FirstTime( PrevSimTime, PrevClockTime, SimStrtTime, SimStrtCPU, t, tMax ) + CALL SimStatus_FirstTime( PrevSimTime, PrevClockTime, SimStrtTime, SimStrtCPU, t, TMax ) DO i = 1,nt @@ -563,7 +566,7 @@ PROGRAM MoorDyn_Driver if ( MOD( i, 20 ) == 0 ) THEN - CALL SimStatus( PrevSimTime, PrevClockTime, t, tMax ) + CALL SimStatus( PrevSimTime, PrevClockTime, t, TMax ) end if ! shift older inputs back in the buffer @@ -572,16 +575,20 @@ PROGRAM MoorDyn_Driver MD_uTimes(2) = MD_uTimes(1) - dtC !MD_uTimes(3) = MD_uTimes(2) - dtC - ! update coupled object kinematics iff we're reading input time series + ! update coupled object kinematics if we're reading input time series if (drvrInitInp%InputsMod == 1 ) then DO iTurb = 1, MD_p%nTurbines K = 1 ! the index of the coupling points in the input mesh CoupledKinematics J = 1 ! the starting index of the relevant DOFs in the input array + IF (MD_InitInp%FarmSize < 0) THEN + MD_p%TurbineRefPos(:,iTurb) = 0.0 + ENDIF + ! any coupled bodies (type -1) DO l = 1,MD_p%nCpldBodies(iTurb) - MD_u(1)%CoupledKinematics(iTurb)%TranslationDisp(:,K) = r_in(i, J:J+2) - MD_u(1)%CoupledKinematics(iTurb)%Position(:,K) - MD_p%TurbineRefPos(:,iTurb) + MD_u(1)%CoupledKinematics(iTurb)%TranslationDisp(:,K) = r_in(i, J:J+2) - MD_u(1)%CoupledKinematics(iTurb)%Position(:,K) - MD_p%TurbineRefPos(:,iTurb) MD_u(1)%CoupledKinematics(iTurb)%Orientation( :,:,K) = EulerConstruct( r_in(i, J+3:J+5) ) ! full Euler angle approach MD_u(1)%CoupledKinematics(iTurb)%TranslationVel( :,K) = rd_in(i, J:J+2) MD_u(1)%CoupledKinematics(iTurb)%RotationVel( :,K) = rd_in(i, J+3:J+5) @@ -594,8 +601,7 @@ PROGRAM MoorDyn_Driver ! any coupled rods (type -1 or -2) >>> need to make rotations ignored if it's a pinned rod <<< DO l = 1,MD_p%nCpldRods(iTurb) - - MD_u(1)%CoupledKinematics(iTurb)%TranslationDisp(:,K) = r_in(i, J:J+2) - MD_u(1)%CoupledKinematics(iTurb)%Position(:,K) - MD_p%TurbineRefPos(:,iTurb) + MD_u(1)%CoupledKinematics(iTurb)%TranslationDisp(:,K) = r_in(i, J:J+2) - MD_u(1)%CoupledKinematics(iTurb)%Position(:,K) - MD_p%TurbineRefPos(:,iTurb) MD_u(1)%CoupledKinematics(iTurb)%Orientation( :,:,K) = EulerConstruct( r_in(i, J+3:J+5) ) MD_u(1)%CoupledKinematics(iTurb)%TranslationVel( :,K) = rd_in(i, J:J+2) MD_u(1)%CoupledKinematics(iTurb)%RotationVel( :,K) = rd_in(i, J+3:J+5) @@ -608,8 +614,7 @@ PROGRAM MoorDyn_Driver ! any coupled points (type -1) DO l = 1, MD_p%nCpldPoints(iTurb) - - MD_u(1)%CoupledKinematics(iTurb)%TranslationDisp(:,K) = r_in(i, J:J+2) - MD_u(1)%CoupledKinematics(iTurb)%Position(:,K) - MD_p%TurbineRefPos(:,iTurb) + MD_u(1)%CoupledKinematics(iTurb)%TranslationDisp(:,K) = r_in(i, J:J+2) - MD_u(1)%CoupledKinematics(iTurb)%Position(:,K) - MD_p%TurbineRefPos(:,iTurb) MD_u(1)%CoupledKinematics(iTurb)%TranslationVel( :,K) = rd_in(i, J:J+2) MD_u(1)%CoupledKinematics(iTurb)%TranslationAcc( :,K) = 0.0_DbKi !rdd_in(i, J:J+2) @@ -669,6 +674,11 @@ PROGRAM MoorDyn_Driver call MD_DestroyInput( MD_u(j), ErrStat2, ErrMsg2) end do + if ( ErrStat /= ErrID_None ) THEN ! Display all errors + CALL WrScr1( "Errors: " ) + CALL WrScr( trim(GetErrStr(ErrStat))//': '//trim(ErrMsg) ) + endif + !close (un) call CleanUp() CALL NormStop() @@ -683,8 +693,8 @@ SUBROUTINE AbortIfFailed() if (ErrStat >= AbortErrLev) then call CleanUp() Call ProgAbort(trim(ErrMsg)) - elseif ( ErrStat /= ErrID_None ) THEN - CALL WrScr1( trim(GetErrStr(ErrStat))//': '//trim(ErrMsg) ) + elseif ( ErrStat2 /= ErrID_None ) THEN + CALL WrScr1( trim(GetErrStr(ErrStat2))//': '//trim(ErrMsg2)//NewLine) end if END SUBROUTINE AbortIfFailed @@ -710,7 +720,7 @@ SUBROUTINE ReadDriverInputFile( inputFile, InitInp) ! Local variables INTEGER :: J ! generic integer for counting - CHARACTER(1024) :: EchoFile ! Name of MoorDyn echo file + ! CHARACTER(1024) :: EchoFile ! Name of MoorDyn echo file CHARACTER(1024) :: FileName ! Name of MoorDyn input file CHARACTER(1024) :: FilePath ! Name of path to MoorDyn input file @@ -728,17 +738,17 @@ SUBROUTINE ReadDriverInputFile( inputFile, InitInp) ! Read until "echo" CALL ReadCom( UnIn, FileName, 'MoorDyn Driver input file header line 1', ErrStat2, ErrMsg2); call AbortIfFailed() CALL ReadCom( UnIn, FileName, 'MoorDyn Driver input file header line 2', ErrStat2, ErrMsg2); call AbortIfFailed() - CALL ReadVar ( UnIn, FileName, InitInp%Echo, 'Echo', 'Echo Input', ErrStat2, ErrMsg2); call AbortIfFailed() - ! If we echo, we rewind - IF ( InitInp%Echo ) THEN - EchoFile = TRIM(FileName)//'.echo' - CALL GetNewUnit( UnEcho ) - CALL OpenEcho ( UnEcho, EchoFile, ErrStat2, ErrMsg2 ); call AbortIfFailed() - REWIND(UnIn) - CALL ReadCom( UnIn, FileName, 'MoorDyn Driver input file header line 1', ErrStat2, ErrMsg2, UnEcho); call AbortIfFailed() - CALL ReadCom( UnIn, FileName, 'MoorDyn Driver input file header line 2', ErrStat2, ErrMsg2, UnEcho); call AbortIfFailed() - CALL ReadVar ( UnIn, FileName, InitInp%Echo, 'Echo', 'Echo the input file data', ErrStat2, ErrMsg2, UnEcho); call AbortIfFailed() - END IF + ! CALL ReadVar ( UnIn, FileName, InitInp%Echo, 'Echo', 'Echo Input', ErrStat2, ErrMsg2); call AbortIfFailed() + ! ! If we echo, we rewind + ! IF ( InitInp%Echo ) THEN + ! EchoFile = TRIM(FileName)//'.echo' + ! CALL GetNewUnit( UnEcho ) + ! CALL OpenEcho ( UnEcho, EchoFile, ErrStat2, ErrMsg2 ); call AbortIfFailed() + ! REWIND(UnIn) + ! CALL ReadCom( UnIn, FileName, 'MoorDyn Driver input file header line 1', ErrStat2, ErrMsg2, UnEcho); call AbortIfFailed() + ! CALL ReadCom( UnIn, FileName, 'MoorDyn Driver input file header line 2', ErrStat2, ErrMsg2, UnEcho); call AbortIfFailed() + ! CALL ReadVar ( UnIn, FileName, InitInp%Echo, 'Echo', 'Echo the input file data', ErrStat2, ErrMsg2, UnEcho); call AbortIfFailed() + ! END IF !---------------------- ENVIRONMENTAL CONDITIONS ------------------------------------------------- CALL ReadCom( UnIn, FileName, 'Environmental conditions header', ErrStat2, ErrMsg2, UnEcho); call AbortIfFailed() CALL ReadVar( UnIn, FileName, InitInp%Gravity, 'Gravity', 'Gravity', ErrStat2, ErrMsg2, UnEcho); call AbortIfFailed() diff --git a/modules/moordyn/src/MoorDyn_IO.f90 b/modules/moordyn/src/MoorDyn_IO.f90 index 271d9791d9..245c9afa1f 100644 --- a/modules/moordyn/src/MoorDyn_IO.f90 +++ b/modules/moordyn/src/MoorDyn_IO.f90 @@ -242,8 +242,7 @@ SUBROUTINE getCoefficientOrCurve(inputString, LineProp_c, LineProp_npoints, Line LineProp_npoints = 0; else ! otherwise interpet the input as a file name to load stress-strain lookup data from - - CALL WrScr("found A letter in the line coefficient value so will try to load the filename.") + CALL WrScr1(" Found a letter in the line EA coefficient value so will try to load the filename.") LineProp_c = 0.0 @@ -251,8 +250,13 @@ SUBROUTINE getCoefficientOrCurve(inputString, LineProp_c, LineProp_npoints, Line CALL GetNewUnit( UnCoef ) CALL OpenFInpFile( UnCoef, TRIM(inputString), ErrStat4, ErrMsg4 ) ! add error handling? + IF (ErrStat4 == ErrID_Fatal) then + ErrStat3 = ErrStat4 + ErrMsg3 = ErrMsg4 + RETURN + ENDIF - READ(UnCoef,'(A)',IOSTAT=ErrStat4) Line2 ! skip the first two lines (title, names, and units) then parse + READ(UnCoef,'(A)',IOSTAT=ErrStat4) Line2 ! skip the first three lines (title, names, and units) then parse READ(UnCoef,'(A)',IOSTAT=ErrStat4) Line2 READ(UnCoef,'(A)',IOSTAT=ErrStat4) Line2 diff --git a/modules/moordyn/src/MoorDyn_Line.f90 b/modules/moordyn/src/MoorDyn_Line.f90 index 37f406486a..190cc4d7eb 100644 --- a/modules/moordyn/src/MoorDyn_Line.f90 +++ b/modules/moordyn/src/MoorDyn_Line.f90 @@ -213,25 +213,6 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) RETURN END IF - - if (p%writeLog > 1) then - write(p%UnLog, '(A)') " - Line"//trim(num2lstr(Line%IdNum)) - write(p%UnLog, '(A)') " ID: "//trim(num2lstr(Line%IdNum)) - write(p%UnLog, '(A)') " UnstrLen: "//trim(num2lstr(Line%UnstrLen)) - write(p%UnLog, '(A)') " N : "//trim(num2lstr(Line%N )) - write(p%UnLog, '(A)') " d : "//trim(num2lstr(Line%d )) - write(p%UnLog, '(A)') " rho : "//trim(num2lstr(Line%rho )) - write(p%UnLog, '(A)') " E : "//trim(num2lstr(Line%EA )) - write(p%UnLog, '(A)') " EI : "//trim(num2lstr(Line%EI )) - !write(p%UnLog, '(A)') " BAin: "//trim(num2lstr(Line%BAin)) - write(p%UnLog, '(A)') " Can : "//trim(num2lstr(Line%Can )) - write(p%UnLog, '(A)') " Cat : "//trim(num2lstr(Line%Cat )) - write(p%UnLog, '(A)') " Cdn : "//trim(num2lstr(Line%Cdn )) - write(p%UnLog, '(A)') " Cdt : "//trim(num2lstr(Line%Cdt )) - !write(p%UnLog, '(A)') " ww_l: " << ( (rho - env->rho_w)*(pi/4.*d*d) )*9.81 << endl; - end if - - ! need to add cleanup sub <<< diff --git a/modules/moordyn/src/MoorDyn_Misc.f90 b/modules/moordyn/src/MoorDyn_Misc.f90 index ddc8bf25a6..8fb91cf218 100644 --- a/modules/moordyn/src/MoorDyn_Misc.f90 +++ b/modules/moordyn/src/MoorDyn_Misc.f90 @@ -151,8 +151,10 @@ subroutine GetOrientationAngles(vec, phi, sinPhi, cosPhi, tanPhi, beta, sinBeta, vecLen = SQRT(Dot_Product(vec,vec)) vecLen2D = SQRT(vec(1)**2+vec(2)**2) if ( vecLen < 0.000001 ) then - print *, "ERROR in GetOrientationAngles in MoorDyn. Supplied vector is near zero" - print *, vec + if (wordy > 0) then + print *, "ERROR in GetOrientationAngles in MoorDyn. Supplied vector is near zero" + print *, vec + endif k_hat = NaN ! 1.0/0.0 else k_hat = vec / vecLen @@ -1348,7 +1350,7 @@ SUBROUTINE setupWaterKin(WaterKinString, p, Tmax, ErrStat, ErrMsg) ELSE IF (SCAN(WaterKinString, "abcdfghijklmnopqrstuvwxyzABCDFGHIJKLMNOPQRSTUVWXYZ") == 0) THEN ! If the input has no letters, let's assume it's a number - print *, "ERROR WaveKin option does not currently support numeric entries. It must be a filename." + call WrScr( "ERROR WaveKin option does not currently support numeric entries. It must be a filename." ) p%WaveKin = 0 p%Current = 0 return @@ -1356,7 +1358,7 @@ SUBROUTINE setupWaterKin(WaterKinString, p, Tmax, ErrStat, ErrMsg) ! otherwise interpret the input as a file name to load the bathymetry lookup data from - print *, " The waterKin input contains letters so will load a water kinematics input file" + call WrScr( " The waterKin input contains letters so will load a water kinematics input file" ) ! -------- load water kinematics input file ------------- @@ -1387,14 +1389,17 @@ SUBROUTINE setupWaterKin(WaterKinString, p, Tmax, ErrStat, ErrMsg) READ(UnIn,*, IOSTAT=ErrStat2) coordtype ! get the entry type READ(UnIn,'(A)', IOSTAT=ErrStat2) entries2 ! get entries as string to be processed CALL gridAxisCoords(coordtype, entries2, p%pxWave, p%nxWave, ErrStat2, ErrMsg2) + Call SetErrStat(ErrStat2,ErrMsg2, ErrStat, ErrMsg, 'MD_getWaterKin') ! Y grid points READ(UnIn,*, IOSTAT=ErrStat2) coordtype ! get the entry type READ(UnIn,'(A)', IOSTAT=ErrStat2) entries2 ! get entries as string to be processed CALL gridAxisCoords(coordtype, entries2, p%pyWave, p%nyWave, ErrStat2, ErrMsg2) + Call SetErrStat(ErrStat2,ErrMsg2, ErrStat, ErrMsg, 'MD_getWaterKin') ! Z grid points READ(UnIn,*, IOSTAT=ErrStat2) coordtype ! get the entry type READ(UnIn,'(A)', IOSTAT=ErrStat2) entries2 ! get entries as string to be processed CALL gridAxisCoords(coordtype, entries2, p%pzWave, p%nzWave, ErrStat2, ErrMsg2) + Call SetErrStat(ErrStat2,ErrMsg2, ErrStat, ErrMsg, 'MD_getWaterKin') ! ----- current ----- CALL ReadCom( UnIn, FileName, 'current header', ErrStat2, ErrMsg2, UnEcho); if(Failed()) return CALL ReadVar( UnIn, FileName, p%Current, 'CurrentMod', 'CurrentMod', ErrStat2, ErrMsg2, UnEcho); if(Failed()) return @@ -1408,7 +1413,7 @@ SUBROUTINE setupWaterKin(WaterKinString, p, Tmax, ErrStat, ErrMsg) EXIT ! break out of the loop if it couldn't read the line (i.e. if at end of file) end if if (i == 100) then - print*,"WARNING: MD can handle a maximum of 100 current profile points" + call WrScr("WARNING: MD can handle a maximum of 100 current profile points") exit end if END DO @@ -1444,7 +1449,7 @@ SUBROUTINE setupWaterKin(WaterKinString, p, Tmax, ErrStat, ErrMsg) ! --------------------- set from inputted wave elevation time series, grid approach ------------------- if (p%WaveKin == 3) then - print *, 'Setting up WaveKin 3 option: read wave elevation time series from file' + call WrScr( 'Setting up WaveKin 3 option: read wave elevation time series from file' ) IF ( LEN_TRIM( WaveKinFile ) == 0 ) THEN CALL SetErrStat( ErrID_Fatal,'WaveKinFile must not be an empty string.',ErrStat, ErrMsg, RoutineName); return @@ -1462,7 +1467,7 @@ SUBROUTINE setupWaterKin(WaterKinString, p, Tmax, ErrStat, ErrMsg) CALL OpenFInpFile ( UnElev, WaveKinFile, ErrStat2, ErrMsg2 ); if(Failed()) return - print *, 'Reading wave elevation data from ', trim(WaveKinFile) + call WrScr( 'Reading wave elevation data from '//trim(WaveKinFile) ) ! Read through length of file to find its length i = 1 ! start counter @@ -1497,7 +1502,7 @@ SUBROUTINE setupWaterKin(WaterKinString, p, Tmax, ErrStat, ErrMsg) ! Close the inputs file CLOSE ( UnElev ) - print *, "Read ", ntIn, " time steps from input file." + call WrScr( "Read "//trim(num2lstr(ntIn))//" time steps from input file." ) ! if (WaveTimeIn(ntIn) < TMax) then <<<< need to handle if time series is too short? @@ -1707,49 +1712,58 @@ SUBROUTINE gridAxisCoords(coordtype, entries, coordarray, n, ErrStat, ErrMsg) REAL(ReKi) :: tempArray (100) REAL(ReKi) :: dx INTEGER(IntKi) :: nEntries, I - - ! get array of coordinate entries - CALL stringToArray(entries, nEntries, tempArray) - - ! set number of coordinates - if ( coordtype==0) then ! 0: not used - make one grid point at zero - n = 1; - else if (coordtype==1) then ! 1: list values in ascending order - n = nEntries - else if (coordtype==2) then ! 2: uniform specified by -xlim, xlim, num - n = int(tempArray(3)) - else - print *, "Error: invalid coordinate type specified to gridAxisCoords" - end if - - ! allocate coordinate array - CALL AllocAry(coordarray, n, 'x,y, or z grid points' , ErrStat, ErrMsg) - !ALLOCATE ( coordarray(n), STAT=ErrStat) - - ! fill in coordinates - if ( coordtype==0) then - coordarray(1) = 0.0_ReKi - - else if (coordtype==1) then - coordarray(1:n) = tempArray(1:n) - - else if (coordtype==2) then - coordarray(1) = tempArray(1) - coordarray(n) = tempArray(2) - dx = (coordarray(n)-coordarray(0))/REAL(n-1) - do i=2,n-1 - coordarray(i) = coordarray(1) + REAL(i)*dx - end do - - else - print *, "Error: invalid coordinate type specified to gridAxisCoords" - end if - - print *, "Set water grid coordinates to :" - DO i=1,n - print *, " ", coordarray(i) - end do - + + IF (len(trim(entries)) == len(entries)) THEN + call WrScr("Warning: Only 120 characters read from wave grid coordinates") + END IF + + IF (entries(len(entries):len(entries)) == ',') THEN + ErrStat = ErrID_Fatal + ErrMsg = 'Last character of wave grid coordinate list cannot be comma' + ELSE + ! get array of coordinate entries + CALL stringToArray(entries, nEntries, tempArray) + + ! set number of coordinates + if ( coordtype==0) then ! 0: not used - make one grid point at zero + n = 1; + else if (coordtype==1) then ! 1: list values in ascending order + n = nEntries + else if (coordtype==2) then ! 2: uniform specified by -xlim, xlim, num + n = int(tempArray(3)) + else + call WrScr("Error: invalid coordinate type specified to gridAxisCoords") + end if + + ! allocate coordinate array + CALL AllocAry(coordarray, n, 'x,y, or z grid points' , ErrStat, ErrMsg) + !ALLOCATE ( coordarray(n), STAT=ErrStat) + + ! fill in coordinates + if ( coordtype==0) then + coordarray(1) = 0.0_ReKi + + else if (coordtype==1) then + coordarray(1:n) = tempArray(1:n) + + else if (coordtype==2) then + coordarray(1) = tempArray(1) + coordarray(n) = tempArray(2) + dx = (coordarray(n)-coordarray(1))/REAL(n-1) + do i=2,n + coordarray(i) = coordarray(i-1) + dx + end do + + else + call WrScr("Error: invalid coordinate type specified to gridAxisCoords") + end if + + ! print *, "Set water grid coordinates to :" + ! DO i=1,n + ! print *, " ", coordarray(i) + ! end do + END IF + END SUBROUTINE gridAxisCoords @@ -1777,7 +1791,7 @@ SUBROUTINE stringToArray(instring, n, outarray) END IF n = n + 1 if (n > 100) then - print *, "ERROR - stringToArray cannot do more than 100 entries" + call WrScr( "ERROR - stringToArray cannot do more than 100 entries") end if READ(instring(pos1:pos1+pos2-2), *) outarray(n) diff --git a/modules/moordyn/src/MoorDyn_Point.f90 b/modules/moordyn/src/MoorDyn_Point.f90 index f79d86858b..9328f4805d 100644 --- a/modules/moordyn/src/MoorDyn_Point.f90 +++ b/modules/moordyn/src/MoorDyn_Point.f90 @@ -362,7 +362,7 @@ SUBROUTINE Point_AddLine(Point, lineID, TopOfLine) Point%Attached(Point%nAttached) = lineID Point%Top(Point%nAttached) = TopOfLine ! attached to line ... 1 = top/fairlead(end B), 0 = bottom/anchor(end A) ELSE - Print*, "Too many lines connected to Point ", Point%IdNum, " in MoorDyn!" + call WrScr("Too many lines connected to Point "//trim(num2lstr(Point%IdNum))//" in MoorDyn!") END IF END SUBROUTINE Point_AddLine @@ -399,7 +399,7 @@ SUBROUTINE Point_RemoveLine(Point, lineID, TopOfLine, rEnd, rdEnd) rdEnd(J) = Point%rd(J) END DO - print*, "Detached line ", lineID, " from Point ", Point%IdNum + call WrScr( "Detached line "//trim(num2lstr(lineID))//" from Point "//trim(num2lstr(Point%IdNum))) EXIT END DO diff --git a/modules/moordyn/src/MoorDyn_Registry.txt b/modules/moordyn/src/MoorDyn_Registry.txt index afe0613139..d4df4982c8 100644 --- a/modules/moordyn/src/MoorDyn_Registry.txt +++ b/modules/moordyn/src/MoorDyn_Registry.txt @@ -24,7 +24,7 @@ typedef MoorDyn/MD InitInputType ReKi g - -99 typedef ^ ^ ReKi rhoW - -999.9 - "sea density" "[kg/m^3]" typedef ^ ^ ReKi WtrDepth - -999.9 - "depth of water" "[m]" typedef ^ ^ ReKi PtfmInit {:}{:} - - "initial position of platform(s) shape: 6, nTurbines" - -typedef ^ ^ IntKi FarmSize - 0 - "Indicates normal FAST module mode if 0, FAST.Farm coupled mode and =nTurbines if >0" - +typedef ^ ^ IntKi FarmSize - 0 - "Indicates normal FAST module mode if 0, FAST.Farm coupled mode and =nTurbines if >0, standalone mode if -1" - typedef ^ ^ ReKi TurbineRefPos {:}{:} - - "reference position of turbines in farm, shape: 3, nTurbines" - typedef ^ ^ ReKi Tmax - - - "simulation duration" "[s]" typedef ^ ^ CHARACTER(1024) FileName - "" - "MoorDyn input file" @@ -87,7 +87,7 @@ typedef ^ ^ DbKi CaEnd - # this is the Body type, which holds data for each body object typedef ^ MD_Body IntKi IdNum - - - "integer identifier of this Point" -typedef ^ ^ IntKi typeNum - - - "integer identifying the type. 0=free, 1=fixed, -1=vessel" +typedef ^ ^ IntKi typeNum - - - "integer identifying the type. 0=free, 1=fixed, -1=coupled, 2=coupledpinned" typedef ^ ^ IntKi AttachedC {30} - - "list of IdNums of points attached to this body" typedef ^ ^ IntKi AttachedR {30} - - "list of IdNums of rods attached to this body" typedef ^ ^ IntKi nAttachedC - 0 - "number of attached points" @@ -116,7 +116,7 @@ typedef ^ ^ DbKi rCG {3} # this is the Point type, which holds data for each point object typedef ^ MD_Point IntKi IdNum - - - "integer identifier of this point" typedef ^ ^ CHARACTER(10) type - - - "type of point: fix, vessel, point" -typedef ^ ^ IntKi typeNum - - - "integer identifying the type. 1=fixed, -1=vessel, 0=free" +typedef ^ ^ IntKi typeNum - - - "integer identifying the type. 1=fixed, -1=coupled, 0=free" typedef ^ ^ IntKi Attached {10} - - "list of IdNums of lines attached to this point node" typedef ^ ^ IntKi Top {10} - - "list of ints specifying whether each line is attached at 1 = top/fairlead(end B), 0 = bottom/anchor(end A)" typedef ^ ^ IntKi nAttached - 0 - "number of attached lines" @@ -142,7 +142,7 @@ typedef ^ ^ DbKi M {3}{3} typedef ^ MD_Rod IntKi IdNum - - - "integer identifier of this Line" typedef ^ ^ CHARACTER(10) type - - - "type of Rod. should match one of RodProp names" typedef ^ ^ IntKi PropsIdNum - - - "the IdNum of the associated rod properties" - -typedef ^ ^ IntKi typeNum - - - "integer identifying the type. 0=fixed, 1=vessel, 2=point" +typedef ^ ^ IntKi typeNum - - - "integer identifying the type. 0=free, 1=pinned, 2=fixed, -1=coupledpinned, -2=coupled" typedef ^ ^ IntKi AttachedA {10} - - "list of IdNums of lines attached to end A" typedef ^ ^ IntKi AttachedB {10} - - "list of IdNums of lines attached to end B" typedef ^ ^ IntKi TopA {10} - - "list of ints specifying whether each line is attached at 1 = top/fairlead(end B), 0 = bottom/anchor(end A)" @@ -390,6 +390,8 @@ typedef ^ ^ DbKi mu_kT - typedef ^ ^ DbKi mu_kA - - - "axial kinetic friction coefficient" "(-)" typedef ^ ^ DbKi mc - - - "ratio of the static friction coefficient to the kinetic friction coefficient" "(-)" typedef ^ ^ DbKi cv - - - "saturated damping coefficient" "(-)" +typedef ^ ^ IntKi Standalone - - - "Indicates MoorDyn run as standalone code if 1, coupled if 0" - +typedef ^ ^ IntKi inertialF - 0 - "Indicates MoorDyn returning inertial moments for coupled 6DOF objects. 1 if yes, 0 if no" - # --- parameters for wave and current --- typedef ^ ^ IntKi nxWave - - - "number of x wave grid points" - typedef ^ ^ IntKi nyWave - - - "number of y wave grid points" - diff --git a/modules/moordyn/src/MoorDyn_Rod.f90 b/modules/moordyn/src/MoorDyn_Rod.f90 index f186891ebb..f5f718198a 100644 --- a/modules/moordyn/src/MoorDyn_Rod.f90 +++ b/modules/moordyn/src/MoorDyn_Rod.f90 @@ -151,24 +151,6 @@ SUBROUTINE Rod_Setup(Rod, RodProp, endCoords, p, ErrStat, ErrMsg) IF (wordy > 0) print *, "Set up Rod ",Rod%IdNum, ", type ", Rod%typeNum - - if (p%writeLog > 1) then - write(p%UnLog, '(A)') " - Rod "//trim(num2lstr(Rod%IdNum)) - write(p%UnLog, '(A)') " ID: "//trim(num2lstr(Rod%IdNum)) - write(p%UnLog, '(A)') " UnstrLen: "//trim(num2lstr(Rod%UnstrLen)) - write(p%UnLog, '(A)') " N : "//trim(num2lstr(Rod%N )) - write(p%UnLog, '(A)') " d : "//trim(num2lstr(Rod%d )) - write(p%UnLog, '(A)') " rho : "//trim(num2lstr(Rod%rho )) - write(p%UnLog, '(A)') " Can : "//trim(num2lstr(Rod%Can )) - write(p%UnLog, '(A)') " Cat : "//trim(num2lstr(Rod%Cat )) - write(p%UnLog, '(A)') " CaEnd: "//trim(num2lstr(Rod%CaEnd )) - write(p%UnLog, '(A)') " Cdn : "//trim(num2lstr(Rod%Cdn )) - write(p%UnLog, '(A)') " Cdt : "//trim(num2lstr(Rod%Cdt )) - write(p%UnLog, '(A)') " CdEnd: "//trim(num2lstr(Rod%CdEnd )) - !write(p%UnLog, '(A)') " ww_l: " << ( (rho - env->rho_w)*(pi/4.*d*d) )*9.81 << endl; - end if - - ! need to add cleanup sub <<< @@ -491,8 +473,7 @@ SUBROUTINE Rod_GetStateDeriv(Rod, Xd, m, p) ELSE ! pinned rod, 6 states (rotational only) ! account for moment in response to end A acceleration due to inertial coupling (off-diagonal sub-matrix terms) - !Fnet(4:6) = Fnet(4:6) - MATMUL(M_out(4:6,1:3), Rod%a6(1:3)) ! <<