From fac429a666c1b59d5f0223a07a01ac5d77fe3122 Mon Sep 17 00:00:00 2001 From: Ganesh Vijayakumar Date: Tue, 9 May 2017 18:48:46 -0600 Subject: [PATCH 1/8] BUGFIX: Fixing various bugs related to the setting of velocity at the force nodes and interpolation to the Aerodyn nodes --- glue-codes/fast-cpp/src/OpenFAST.H | 2 +- glue-codes/fast-cpp/src/OpenFAST.cpp | 64 ++++++++++++++++++------- modules-local/openfoam/src/OpenFOAM.f90 | 53 ++++++++++---------- 3 files changed, 77 insertions(+), 42 deletions(-) diff --git a/glue-codes/fast-cpp/src/OpenFAST.H b/glue-codes/fast-cpp/src/OpenFAST.H index aa43032718..709fb0bee7 100644 --- a/glue-codes/fast-cpp/src/OpenFAST.H +++ b/glue-codes/fast-cpp/src/OpenFAST.H @@ -1,6 +1,5 @@ #ifndef OpenFAST_h #define OpenFAST_h - #include "FAST_Library.h" #include "sys/stat.h" #include "math.h" @@ -152,6 +151,7 @@ class OpenFAST { void setTurbineProcNo(int iTurbGlob, int procNo) { turbineMapGlobToProc[iTurbGlob] = procNo; } void allocateTurbinesToProcsSimple(); + void getApproxHubPos(std::vector & currentCoords, int iTurbGlob); void getHubPos(std::vector & currentCoords, int iTurbGlob); void getHubShftDir(std::vector & hubShftVec, int iTurbGlob); diff --git a/glue-codes/fast-cpp/src/OpenFAST.cpp b/glue-codes/fast-cpp/src/OpenFAST.cpp index 5e313d65b1..7cc53b0a43 100644 --- a/glue-codes/fast-cpp/src/OpenFAST.cpp +++ b/glue-codes/fast-cpp/src/OpenFAST.cpp @@ -86,9 +86,9 @@ void fast::OpenFAST::solution0() { if (!dryRun) { // set wind speeds at initial locations - for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { - setOutputsToFAST(cDriver_Input_from_FAST[iTurb], cDriver_Output_to_FAST[iTurb]); - } + // for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { + // setOutputsToFAST(cDriver_Input_from_FAST[iTurb], cDriver_Output_to_FAST[iTurb]); + // } if(scStatus) { @@ -128,7 +128,7 @@ void fast::OpenFAST::step() { for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { // set wind speeds at original locations - setOutputsToFAST(cDriver_Input_from_FAST[iTurb], cDriver_Output_to_FAST[iTurb]); + // setOutputsToFAST(cDriver_Input_from_FAST[iTurb], cDriver_Output_to_FAST[iTurb]); // this advances the states, calls CalcOutput, and solves for next inputs. Predictor-corrector loop is imbeded here: // (note OpenFOAM could do subcycling around this step) @@ -138,15 +138,27 @@ void fast::OpenFAST::step() { } } + if ( isDebug() ) { + + std::ofstream fastcpp_velocity_file; + fastcpp_velocity_file.open("fastcpp_velocity.csv") ; + fastcpp_velocity_file << "# x, y, z, fx, fy, fz" << std::endl ; + fastcpp_velocity_file << get_numVelPtsLoc(iTurb) << std::endl ; + for (int iNode=0; iNode < get_numVelPtsLoc(iTurb); iNode++) { + fastcpp_velocity_file << cDriver_Input_from_FAST[iTurb].pxForce[iNode] << ", " << cDriver_Input_from_FAST[iTurb].pyForce[iNode] << ", " << cDriver_Input_from_FAST[iTurb].pzForce[iNode] << ", " << cDriver_Output_to_FAST[iTurb].u[iNode] << ", " << cDriver_Output_to_FAST[iTurb].v[iNode] << ", " << cDriver_Output_to_FAST[iTurb].w[iNode] << " " << std::endl ; + } + fastcpp_velocity_file.close() ; + + } + + + FAST_OpFM_Step(&iTurb, &ErrStat, ErrMsg); checkError(ErrStat, ErrMsg); if ( isDebug() ) { - for (int iNode=0; iNode < get_numForcePtsLoc(iTurb); iNode++) { - std::cout << "Node " << iNode << " Position = " << cDriver_Input_from_FAST[iTurb].pxForce[iNode] << " " << cDriver_Input_from_FAST[iTurb].pyForce[iNode] << " " << cDriver_Input_from_FAST[iTurb].pzForce[iNode] << " " << std::endl ; - } std::ofstream actuatorForcesFile; - actuatorForcesFile.open("actuatorForces.csv") ; + actuatorForcesFile.open("actuator_forces.csv") ; actuatorForcesFile << "# x, y, z, fx, fy, fz" << std::endl ; for (int iNode=0; iNode < get_numForcePtsLoc(iTurb); iNode++) { actuatorForcesFile << cDriver_Input_from_FAST[iTurb].pxForce[iNode] << ", " << cDriver_Input_from_FAST[iTurb].pyForce[iNode] << ", " << cDriver_Input_from_FAST[iTurb].pzForce[iNode] << ", " << cDriver_Input_from_FAST[iTurb].fx[iNode] << ", " << cDriver_Input_from_FAST[iTurb].fy[iNode] << ", " << cDriver_Input_from_FAST[iTurb].fz[iNode] << " " << std::endl ; @@ -246,12 +258,22 @@ void fast::OpenFAST::setOutputsToFAST(OpFM_InputType_t cDriver_Input_from_FAST, } +void fast::OpenFAST::getApproxHubPos(std::vector & currentCoords, int iTurbGlob) { + + // Get hub position of Turbine 'iTurbGlob' + currentCoords[0] = globTurbineData[iTurbGlob].TurbineHubPos[0]; + currentCoords[1] = globTurbineData[iTurbGlob].TurbineHubPos[1]; + currentCoords[2] = globTurbineData[iTurbGlob].TurbineHubPos[2]; + +} + void fast::OpenFAST::getHubPos(std::vector & currentCoords, int iTurbGlob) { // Get hub position of Turbine 'iTurbGlob' - currentCoords[0] = globTurbineData[iTurbGlob].TurbineHubPos[0] ; - currentCoords[1] = globTurbineData[iTurbGlob].TurbineHubPos[1] ; - currentCoords[2] = globTurbineData[iTurbGlob].TurbineHubPos[2] ; + int iTurbLoc = get_localTurbNo(iTurbGlob); + currentCoords[0] = cDriver_Input_from_FAST[iTurbLoc].pxVel[0] ; + currentCoords[1] = cDriver_Input_from_FAST[iTurbLoc].pyVel[0] ; + currentCoords[2] = cDriver_Input_from_FAST[iTurbLoc].pzVel[0] ; } @@ -323,10 +345,9 @@ void fast::OpenFAST::setVelocity(std::vector & currentVelocity, int iNod // Set velocity at current node of current turbine - int iTurbLoc = get_localTurbNo(iTurbGlob); for(int j=0; j < iTurbLoc; j++) iNode = iNode - get_numVelPtsLoc(iTurbLoc); - cDriver_Output_to_FAST[iTurbLoc].u[iNode] = currentVelocity[0]; - cDriver_Output_to_FAST[iTurbLoc].v[iNode] = currentVelocity[1]; - cDriver_Output_to_FAST[iTurbLoc].w[iNode] = currentVelocity[2]; - + forceNodeVel[iTurbLoc][iNode][0] = currentVelocity[0]; + forceNodeVel[iTurbLoc][iNode][1] = currentVelocity[1]; + forceNodeVel[iTurbLoc][iNode][2] = currentVelocity[2]; } void fast::OpenFAST::interpolateVel_ForceToVelNodes() { @@ -339,6 +360,15 @@ void fast::OpenFAST::interpolateVel_ForceToVelNodes() { cDriver_Output_to_FAST[iTurb].v[0] = forceNodeVel[iTurb][0][1]; cDriver_Output_to_FAST[iTurb].w[0] = forceNodeVel[iTurb][0][2]; + if ( isDebug() ) { + std::ofstream actuatorVelFile; + actuatorVelFile.open("actuator_velocity.csv") ; + for (int iNode=0; iNode < get_numForcePtsLoc(iTurb); iNode++) { + actuatorVelFile << cDriver_Input_from_FAST[iTurb].pxForce[iNode] << ", " << cDriver_Input_from_FAST[iTurb].pyForce[iNode] << ", " << cDriver_Input_from_FAST[iTurb].pzForce[iNode] << ", " << forceNodeVel[iTurb][iNode][0] << ", " << forceNodeVel[iTurb][iNode][1] << ", " << forceNodeVel[iTurb][iNode][2] << " " << std::endl ; + } + actuatorVelFile.close() ; + } + // Do the blades first int nBlades = get_numBladesLoc(iTurb); for(int iBlade=0; iBlade < nBlades; iBlade++) { @@ -366,7 +396,7 @@ void fast::OpenFAST::interpolateVel_ForceToVelNodes() { ); //Find nearest two force nodes int jForceLower = 0; - while ( ((rDistForce[jForceLower]-rDistVel)*(rDistForce[jForceLower+1]-rDistVel) > 0) && (jForceLower < (nForcePtsBlade-1) ) ) { + while ( (rDistForce[jForceLower+1] < rDistVel) && ( jForceLower < (nForcePtsBlade-1)) ) { jForceLower = jForceLower + 1; } int iNodeForceLower = 1 + iBlade * nForcePtsBlade + jForceLower ; @@ -405,7 +435,7 @@ void fast::OpenFAST::interpolateVel_ForceToVelNodes() { ); //Find nearest two force nodes int jForceLower = 0; - while ( ((hDistForce[jForceLower]-hDistVel)*(hDistForce[jForceLower+1]-hDistVel) > 0) && (jForceLower < (nForcePtsTower-1) ) ) { + while ( (hDistForce[jForceLower+1] < hDistVel) && ( jForceLower < (nForcePtsTower-1)) ) { jForceLower = jForceLower + 1; } int iNodeForceLower = iNodeBotTowerForce + jForceLower ; diff --git a/modules-local/openfoam/src/OpenFOAM.f90 b/modules-local/openfoam/src/OpenFOAM.f90 index 6fc0d48caf..a949bb571e 100644 --- a/modules-local/openfoam/src/OpenFOAM.f90 +++ b/modules-local/openfoam/src/OpenFOAM.f90 @@ -476,25 +476,31 @@ SUBROUTINE SetOpFMForces(p_FAST, p_AD14, u_AD14, y_AD14, u_AD, y_AD, y_ED, OpFM, ! blade nodes !....................... - open(unit=987,file='AeroDyn_forces.csv') + open(unit=987,file='fast_aerodyn_velocity_forces.csv') write(987,*) '#x, y, z, fx, fy, fz' + + open(unit=988,file='fast_actuator_forces.csv') + write(988,*) '#x, y, z, fx, fy, fz' + DO K = 1,OpFM%p%NumBl + DO J = 1,u_AD%BladeMotion(k)%NNodes + write(987,*) u_AD%BladeMotion(k)%TranslationDisp(1,j) + u_AD%BladeMotion(k)%Position(1,j), ', ', u_AD%BladeMotion(k)%TranslationDisp(2,j) + u_AD%BladeMotion(k)%Position(2,j), ', ', u_AD%BladeMotion(k)%TranslationDisp(3,j) + u_AD%BladeMotion(k)%Position(3,j), ', ', OpFM%y%u(1 + (k-1)*u_AD%BladeMotion(k)%NNodes + j), ', ', OpFM%y%v(1 + (k-1)*u_AD%BladeMotion(k)%NNodes + j), ', ', OpFM%y%w(1 + (k-1)*u_AD%BladeMotion(k)%NNodes + j), ', ', y_AD%BladeLoad(k)%Force(1,j), ', ', y_AD%BladeLoad(k)%Force(2,j), ', ', y_AD%BladeLoad(k)%Force(2,j) + END DO + call Transfer_Line2_to_Line2( y_AD%BladeLoad(k), OpFM%m%ActForceLoads(k), OpFM%m%Line2_to_Line2_Loads(k), ErrStat2, ErrMsg2, u_AD%BladeMotion(k), OpFM%m%ActForceMotions(k) ) call Transfer_Line2_to_Point( OpFM%m%ActForceLoads(k), OpFM%m%ActForceLoadsPoints(k), OpFM%m%Line2_to_Point_Loads(k), ErrStat2, ErrMsg2, OpFM%m%ActForceMotions(k), OpFM%m%ActForceMotionsPoints(k) ) - DO J = 1,u_AD%BladeMotion(k)%NNodes - write(987,*) u_AD%BladeMotion(k)%TranslationDisp(1,j) + u_AD%BladeMotion(k)%Position(1,j), ', ', u_AD%BladeMotion(k)%TranslationDisp(2,j) + u_AD%BladeMotion(k)%Position(2,j), ', ', u_AD%BladeMotion(k)%TranslationDisp(3,j) + u_AD%BladeMotion(k)%Position(3,j), ', ', y_AD%BladeLoad(k)%Force(1,j), ', ', y_AD%BladeLoad(k)%Force(2,j), ', ', y_AD%BladeLoad(k)%Force(2,j) - END DO - DO J = 1, OpFM%p%NnodesForceBlade Node = Node + 1 - OpFM%u%fx(Node) = OpFM%m%ActForceLoadsPoints(k)%Force(1,j) / OpFM%p%AirDens - OpFM%u%fy(Node) = OpFM%m%ActForceLoadsPoints(k)%Force(2,j) / OpFM%p%AirDens - OpFM%u%fz(Node) = OpFM%m%ActForceLoadsPoints(k)%Force(3,j) / OpFM%p%AirDens - OpFM%u%momentx(Node) = OpFM%m%ActForceLoadsPoints(k)%Moment(1,j) / OpFM%p%AirDens - OpFM%u%momenty(Node) = OpFM%m%ActForceLoadsPoints(k)%Moment(2,j) / OpFM%p%AirDens - OpFM%u%momentz(Node) = OpFM%m%ActForceLoadsPoints(k)%Moment(3,j) / OpFM%p%AirDens + OpFM%u%fx(Node) = OpFM%m%ActForceLoadsPoints(k)%Force(1,j) + OpFM%u%fy(Node) = OpFM%m%ActForceLoadsPoints(k)%Force(2,j) + OpFM%u%fz(Node) = OpFM%m%ActForceLoadsPoints(k)%Force(3,j) + OpFM%u%momentx(Node) = OpFM%m%ActForceLoadsPoints(k)%Moment(1,j) + OpFM%u%momenty(Node) = OpFM%m%ActForceLoadsPoints(k)%Moment(2,j) + OpFM%u%momentz(Node) = OpFM%m%ActForceLoadsPoints(k)%Moment(3,j) + + write(988,*) OpFM%u%pxForce(Node), ', ', OpFM%u%pyForce(Node), ', ', OpFM%u%pzForce(Node), ', ', OpFM%u%fx(Node), ', ', OpFM%u%fy(Node), ', ', OpFM%u%fz(Node), ', ' END DO END DO !K = 1,OpFM%p%NumBl @@ -507,26 +513,26 @@ SUBROUTINE SetOpFMForces(p_FAST, p_AD14, u_AD14, y_AD14, u_AD, y_AD, y_ED, OpFM, k = SIZE(u_AD%BladeMotion) + 1 DO J = 1,u_AD%TowerMotion%NNodes - write(987,*) u_AD%TowerMotion%TranslationDisp(1,j) + u_AD%TowerMotion%Position(1,j), ', ', u_AD%TowerMotion%TranslationDisp(2,j) + u_AD%TowerMotion%Position(2,j), ', ', u_AD%TowerMotion%TranslationDisp(3,j) + u_AD%TowerMotion%Position(3,j), ', ', y_AD%TowerLoad%Force(1,j), ', ', y_AD%TowerLoad%Force(2,j), ', ', y_AD%TowerLoad%Force(2,j) + write(987,*) u_AD%TowerMotion%TranslationDisp(1,j) + u_AD%TowerMotion%Position(1,j), ', ', u_AD%TowerMotion%TranslationDisp(2,j) + u_AD%TowerMotion%Position(2,j), ', ', u_AD%TowerMotion%TranslationDisp(3,j) + u_AD%TowerMotion%Position(3,j), ', ', OpFM%y%u(1 + OpFM%p%NumBl*u_AD%BladeMotion(k)%NNodes + j), ', ', OpFM%y%v(1 + OpFM%p%NumBl*u_AD%BladeMotion(k)%NNodes + j), ', ', OpFM%y%w(1 + OpFM%p%NumBl*u_AD%BladeMotion(k)%NNodes + j), ', ', y_AD%TowerLoad%Force(1,j), ', ', y_AD%TowerLoad%Force(2,j), ', ', y_AD%TowerLoad%Force(2,j) END DO call Transfer_Line2_to_Line2( y_AD%TowerLoad, OpFM%m%ActForceLoads(k), OpFM%m%Line2_to_Line2_Loads(k), ErrStat2, ErrMsg2, u_AD%TowerMotion, OpFM%m%ActForceMotions(k) ) call Transfer_Line2_to_Point( OpFM%m%ActForceLoads(k), OpFM%m%ActForceLoadsPoints(k), OpFM%m%Line2_to_Point_Loads(k), ErrStat2, ErrMsg2, OpFM%m%ActForceMotions(k), OpFM%m%ActForceMotionsPoints(k) ) - DO J = 1,u_AD%TowerMotion%NNodes - write(987,*) u_AD%TowerMotion%TranslationDisp(1,j) + u_AD%TowerMotion%Position(1,j), ', ', u_AD%TowerMotion%TranslationDisp(2,j) + u_AD%TowerMotion%Position(2,j), ', ', u_AD%TowerMotion%TranslationDisp(3,j) + u_AD%TowerMotion%Position(3,j), ', ', y_AD%TowerLoad%Force(1,j), ', ', y_AD%TowerLoad%Force(2,j), ', ', y_AD%TowerLoad%Force(2,j) - END DO - DO J=1,OpFM%p%NnodesForceTower Node = Node + 1 - OpFM%u%fx(Node) = OpFM%m%ActForceLoadsPoints(k)%Force(1,j) / OpFM%p%AirDens - OpFM%u%fy(Node) = OpFM%m%ActForceLoadsPoints(k)%Force(2,j) / OpFM%p%AirDens - OpFM%u%fz(Node) = OpFM%m%ActForceLoadsPoints(k)%Force(3,j) / OpFM%p%AirDens - OpFM%u%momentx(Node) = OpFM%m%ActForceLoadsPoints(k)%Moment(1,j) / OpFM%p%AirDens - OpFM%u%momenty(Node) = OpFM%m%ActForceLoadsPoints(k)%Moment(2,j) / OpFM%p%AirDens - OpFM%u%momentz(Node) = OpFM%m%ActForceLoadsPoints(k)%Moment(3,j) / OpFM%p%AirDens + OpFM%u%fx(Node) = OpFM%m%ActForceLoadsPoints(k)%Force(1,j) + OpFM%u%fy(Node) = OpFM%m%ActForceLoadsPoints(k)%Force(2,j) + OpFM%u%fz(Node) = OpFM%m%ActForceLoadsPoints(k)%Force(3,j) + OpFM%u%momentx(Node) = OpFM%m%ActForceLoadsPoints(k)%Moment(1,j) + OpFM%u%momenty(Node) = OpFM%m%ActForceLoadsPoints(k)%Moment(2,j) + OpFM%u%momentz(Node) = OpFM%m%ActForceLoadsPoints(k)%Moment(3,j) + + write(988,*) OpFM%u%pxForce(Node), ', ', OpFM%u%pyForce(Node), ', ', OpFM%u%pzForce(Node), ', ', OpFM%u%fx(Node), ', ', OpFM%u%fy(Node), ', ', OpFM%u%fz(Node), ', ' END DO + close(987) + close(988) END SUBROUTINE SetOpFMForces !---------------------------------------------------------------------------------------------------------------------------------- @@ -535,8 +541,7 @@ SUBROUTINE OpFM_SetWriteOutput( OpFM ) TYPE(OpenFOAM_Data), INTENT(INOUT) :: OpFM ! data for the OpenFOAM integration module - - ! set the hub-height wind speeds + ! set the hub-height wind speeds IF ( ALLOCATED( OpFM%y%WriteOutput ) ) THEN IF ( ASSOCIATED( OpFM%y%u ) ) then OpFM%y%WriteOutput(1) = OpFM%y%u(1) From 69572cf211af767808eb88ebe555f71ffc64da04 Mon Sep 17 00:00:00 2001 From: Ganesh Vijayakumar Date: Mon, 15 May 2017 11:04:43 -0600 Subject: [PATCH 2/8] BUGFIX: Include turbine base location in coordinates returned to driver program --- glue-codes/fast-cpp/src/OpenFAST.cpp | 162 ++++++++++----------------- 1 file changed, 57 insertions(+), 105 deletions(-) diff --git a/glue-codes/fast-cpp/src/OpenFAST.cpp b/glue-codes/fast-cpp/src/OpenFAST.cpp index 7cc53b0a43..602d3ed569 100644 --- a/glue-codes/fast-cpp/src/OpenFAST.cpp +++ b/glue-codes/fast-cpp/src/OpenFAST.cpp @@ -46,6 +46,10 @@ void fast::OpenFAST::init() { FAST_OpFM_Restart(&iTurb, CheckpointFileRoot[iTurb].data(), &AbortErrLev, &dtFAST, &numBlades[iTurb], &numVelPtsBlade[iTurb], &ntStart, &cDriver_Input_from_FAST[iTurb], &cDriver_Output_to_FAST[iTurb], &cDriverSC_Input_from_FAST[iTurb], &cDriverSC_Output_to_FAST[iTurb], &ErrStat, ErrMsg); checkError(ErrStat, ErrMsg); nt_global = ntStart; + + int nfpts = get_numForcePtsLoc(iTurb); + forceNodeVel[iTurb].resize(nfpts); + for (int k = 0; k < nfpts; k++) forceNodeVel[iTurb][k].resize(3) ; } if(scStatus) { @@ -57,7 +61,6 @@ void fast::OpenFAST::init() { // this calls the Init() routines of each module - forceNodeVel.resize(nTurbinesProc); for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { FAST_OpFM_Init(&iTurb, &tMax, FASTInputFileName[iTurb].data(), &TurbID[iTurb], &numScOutputs, &numScInputs, &numForcePtsBlade[iTurb], &numForcePtsTwr[iTurb], TurbineBasePos[iTurb].data(), &AbortErrLev, &dtFAST, &numBlades[iTurb], &numVelPtsBlade[iTurb], &cDriver_Input_from_FAST[iTurb], &cDriver_Output_to_FAST[iTurb], &cDriverSC_Input_from_FAST[iTurb], &cDriverSC_Output_to_FAST[iTurb], &ErrStat, ErrMsg); checkError(ErrStat, ErrMsg); @@ -66,9 +69,9 @@ void fast::OpenFAST::init() { numVelPtsTwr[iTurb] = cDriver_Output_to_FAST[iTurb].u_Len - numBlades[iTurb]*numVelPtsBlade[iTurb] - 1; - int nfPts = get_numForcePtsLoc(iTurb); - forceNodeVel[iTurb].resize(nfPts) ; - for (int k = 0; k < nfPts; k++) forceNodeVel[iTurb][k].resize(3) ; + int nfpts = get_numForcePtsLoc(iTurb); + forceNodeVel[iTurb].resize(nfpts); + for (int k = 0; k < nfpts; k++) forceNodeVel[iTurb][k].resize(3) ; if ( isDebug() ) { for (int iNode=0; iNode < get_numVelPtsLoc(iTurb); iNode++) { @@ -132,26 +135,17 @@ void fast::OpenFAST::step() { // this advances the states, calls CalcOutput, and solves for next inputs. Predictor-corrector loop is imbeded here: // (note OpenFOAM could do subcycling around this step) - if ( isDebug() ) { - for (int iNode=0; iNode < get_numVelPtsLoc(iTurb); iNode++) { - std::cout << "Node " << iNode << " Velocity = " << cDriver_Output_to_FAST[iTurb].u[iNode] << " " << cDriver_Output_to_FAST[iTurb].v[iNode] << " " << cDriver_Output_to_FAST[iTurb].w[iNode] << " " << std::endl ; - } - } - if ( isDebug() ) { std::ofstream fastcpp_velocity_file; fastcpp_velocity_file.open("fastcpp_velocity.csv") ; - fastcpp_velocity_file << "# x, y, z, fx, fy, fz" << std::endl ; - fastcpp_velocity_file << get_numVelPtsLoc(iTurb) << std::endl ; + fastcpp_velocity_file << "# x, y, z, Vx, Vy, Vz" << std::endl ; for (int iNode=0; iNode < get_numVelPtsLoc(iTurb); iNode++) { - fastcpp_velocity_file << cDriver_Input_from_FAST[iTurb].pxForce[iNode] << ", " << cDriver_Input_from_FAST[iTurb].pyForce[iNode] << ", " << cDriver_Input_from_FAST[iTurb].pzForce[iNode] << ", " << cDriver_Output_to_FAST[iTurb].u[iNode] << ", " << cDriver_Output_to_FAST[iTurb].v[iNode] << ", " << cDriver_Output_to_FAST[iTurb].w[iNode] << " " << std::endl ; + fastcpp_velocity_file << cDriver_Input_from_FAST[iTurb].pxVel[iNode] << ", " << cDriver_Input_from_FAST[iTurb].pyVel[iNode] << ", " << cDriver_Input_from_FAST[iTurb].pzVel[iNode] << ", " << cDriver_Output_to_FAST[iTurb].u[iNode] << ", " << cDriver_Output_to_FAST[iTurb].v[iNode] << ", " << cDriver_Output_to_FAST[iTurb].w[iNode] << " " << std::endl ; } fastcpp_velocity_file.close() ; } - - FAST_OpFM_Step(&iTurb, &ErrStat, ErrMsg); checkError(ErrStat, ErrMsg); @@ -271,9 +265,9 @@ void fast::OpenFAST::getHubPos(std::vector & currentCoords, int iTurbGlo // Get hub position of Turbine 'iTurbGlob' int iTurbLoc = get_localTurbNo(iTurbGlob); - currentCoords[0] = cDriver_Input_from_FAST[iTurbLoc].pxVel[0] ; - currentCoords[1] = cDriver_Input_from_FAST[iTurbLoc].pyVel[0] ; - currentCoords[2] = cDriver_Input_from_FAST[iTurbLoc].pzVel[0] ; + currentCoords[0] = cDriver_Input_from_FAST[iTurbLoc].pxVel[0] + TurbineBasePos[iTurbLoc][0] ; + currentCoords[1] = cDriver_Input_from_FAST[iTurbLoc].pyVel[0] + TurbineBasePos[iTurbLoc][1] ; + currentCoords[2] = cDriver_Input_from_FAST[iTurbLoc].pzVel[0] + TurbineBasePos[iTurbLoc][2] ; } @@ -293,9 +287,9 @@ void fast::OpenFAST::getVelNodeCoordinates(std::vector & currentCoords, // Set coordinates at current node of current turbine int iTurbLoc = get_localTurbNo(iTurbGlob); for(int j=0; j < iTurbLoc; j++) iNode = iNode - get_numVelPtsLoc(iTurbLoc); - currentCoords[0] = cDriver_Input_from_FAST[iTurbLoc].pxVel[iNode] ; - currentCoords[1] = cDriver_Input_from_FAST[iTurbLoc].pyVel[iNode] ; - currentCoords[2] = cDriver_Input_from_FAST[iTurbLoc].pzVel[iNode] ; + currentCoords[0] = cDriver_Input_from_FAST[iTurbLoc].pxVel[iNode] + TurbineBasePos[iTurbLoc][0] ; + currentCoords[1] = cDriver_Input_from_FAST[iTurbLoc].pyVel[iNode] + TurbineBasePos[iTurbLoc][1] ; + currentCoords[2] = cDriver_Input_from_FAST[iTurbLoc].pzVel[iNode] + TurbineBasePos[iTurbLoc][2] ; } @@ -303,9 +297,9 @@ void fast::OpenFAST::getForceNodeCoordinates(std::vector & currentCoords // Set coordinates at current node of current turbine int iTurbLoc = get_localTurbNo(iTurbGlob); - currentCoords[0] = cDriver_Input_from_FAST[iTurbLoc].pxForce[iNode] ; - currentCoords[1] = cDriver_Input_from_FAST[iTurbLoc].pyForce[iNode] ; - currentCoords[2] = cDriver_Input_from_FAST[iTurbLoc].pzForce[iNode] ; + currentCoords[0] = cDriver_Input_from_FAST[iTurbLoc].pxForce[iNode] + TurbineBasePos[iTurbLoc][0] ; + currentCoords[1] = cDriver_Input_from_FAST[iTurbLoc].pyForce[iNode] + TurbineBasePos[iTurbLoc][1] ; + currentCoords[2] = cDriver_Input_from_FAST[iTurbLoc].pzForce[iNode] + TurbineBasePos[iTurbLoc][2] ; } @@ -363,6 +357,7 @@ void fast::OpenFAST::interpolateVel_ForceToVelNodes() { if ( isDebug() ) { std::ofstream actuatorVelFile; actuatorVelFile.open("actuator_velocity.csv") ; + actuatorVelFile << "# x, y, z, Vx, Vy, Vz" << std::endl ; for (int iNode=0; iNode < get_numForcePtsLoc(iTurb); iNode++) { actuatorVelFile << cDriver_Input_from_FAST[iTurb].pxForce[iNode] << ", " << cDriver_Input_from_FAST[iTurb].pyForce[iNode] << ", " << cDriver_Input_from_FAST[iTurb].pzForce[iNode] << ", " << forceNodeVel[iTurb][iNode][0] << ", " << forceNodeVel[iTurb][iNode][1] << ", " << forceNodeVel[iTurb][iNode][2] << " " << std::endl ; } @@ -396,7 +391,7 @@ void fast::OpenFAST::interpolateVel_ForceToVelNodes() { ); //Find nearest two force nodes int jForceLower = 0; - while ( (rDistForce[jForceLower+1] < rDistVel) && ( jForceLower < (nForcePtsBlade-1)) ) { + while ( (rDistForce[jForceLower+1] < rDistVel) && ( jForceLower < (nForcePtsBlade-2)) ) { jForceLower = jForceLower + 1; } int iNodeForceLower = 1 + iBlade * nForcePtsBlade + jForceLower ; @@ -435,7 +430,7 @@ void fast::OpenFAST::interpolateVel_ForceToVelNodes() { ); //Find nearest two force nodes int jForceLower = 0; - while ( (hDistForce[jForceLower+1] < hDistVel) && ( jForceLower < (nForcePtsTower-1)) ) { + while ( (hDistForce[jForceLower+1] < hDistVel) && ( jForceLower < (nForcePtsTower-2)) ) { jForceLower = jForceLower + 1; } int iNodeForceLower = iNodeBotTowerForce + jForceLower ; @@ -453,10 +448,14 @@ void fast::OpenFAST::interpolateVel_ForceToVelNodes() { void fast::OpenFAST::computeTorqueThrust(int iTurbGlob, std::vector & torque, std::vector & thrust) { //Compute the torque and thrust based on the forces at the actuator nodes - double relLoc[] = {0.0,0.0,0.0} ; + std::vector relLoc(3,0.0); + std::vector rPerpShft(3); thrust[0] = 0.0; thrust[1] = 0.0; thrust[2] = 0.0; torque[0] = 0.0; torque[1] = 0.0; torque[2] = 0.0; + std::vector hubShftVec(3); + getHubShftDir(hubShftVec, iTurbGlob); + int iTurbLoc = get_localTurbNo(iTurbGlob) ; for (int k=0; k < get_numBladesLoc(iTurbLoc); k++) { for (int j=0; j < numForcePtsBlade[iTurbLoc]; j++) { @@ -469,10 +468,13 @@ void fast::OpenFAST::computeTorqueThrust(int iTurbGlob, std::vector & to relLoc[0] = cDriver_Input_from_FAST[iTurbLoc].pxForce[iNode] - cDriver_Input_from_FAST[iTurbLoc].pxForce[0] ; relLoc[1] = cDriver_Input_from_FAST[iTurbLoc].pyForce[iNode] - cDriver_Input_from_FAST[iTurbLoc].pyForce[0]; relLoc[2] = cDriver_Input_from_FAST[iTurbLoc].pzForce[iNode] - cDriver_Input_from_FAST[iTurbLoc].pzForce[0]; + + double rDotHubShftVec = relLoc[0]*hubShftVec[0] + relLoc[1]*hubShftVec[1] + relLoc[2]*hubShftVec[2]; + for (int j=0; j < 3; j++) rPerpShft[j] = relLoc[j] - rDotHubShftVec * hubShftVec[j]; - torque[0] = torque[0] + relLoc[1] * cDriver_Input_from_FAST[iTurbLoc].fz[iNode] - relLoc[2] * cDriver_Input_from_FAST[iTurbLoc].fy[iNode] + cDriver_Input_from_FAST[iTurbLoc].momentx[iNode] ; - torque[1] = torque[1] + relLoc[2] * cDriver_Input_from_FAST[iTurbLoc].fx[iNode] - relLoc[0] * cDriver_Input_from_FAST[iTurbLoc].fz[iNode] + cDriver_Input_from_FAST[iTurbLoc].momenty[iNode] ; - torque[2] = torque[2] + relLoc[0] * cDriver_Input_from_FAST[iTurbLoc].fy[iNode] - relLoc[1] * cDriver_Input_from_FAST[iTurbLoc].fx[iNode] + cDriver_Input_from_FAST[iTurbLoc].momentz[iNode] ; + torque[0] = torque[0] + rPerpShft[1] * cDriver_Input_from_FAST[iTurbLoc].fz[iNode] - rPerpShft[2] * cDriver_Input_from_FAST[iTurbLoc].fy[iNode] + cDriver_Input_from_FAST[iTurbLoc].momentx[iNode] ; + torque[1] = torque[1] + rPerpShft[2] * cDriver_Input_from_FAST[iTurbLoc].fx[iNode] - rPerpShft[0] * cDriver_Input_from_FAST[iTurbLoc].fz[iNode] + cDriver_Input_from_FAST[iTurbLoc].momenty[iNode] ; + torque[2] = torque[2] + rPerpShft[0] * cDriver_Input_from_FAST[iTurbLoc].fy[iNode] - rPerpShft[1] * cDriver_Input_from_FAST[iTurbLoc].fx[iNode] + cDriver_Input_from_FAST[iTurbLoc].momentz[iNode] ; } } @@ -578,6 +580,7 @@ void fast::OpenFAST::allocateMemory() { numForcePtsTwr.resize(nTurbinesProc); numVelPtsBlade.resize(nTurbinesProc); numVelPtsTwr.resize(nTurbinesProc); + forceNodeVel.resize(nTurbinesProc); for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { @@ -607,57 +610,6 @@ void fast::OpenFAST::allocateMemory() { } -// void fast::OpenFAST::readTurbineData(int iTurb, YAML::Node turbNode) { - -// //Read turbine data for a given turbine using the YAML node -// if (turbNode["turb_id"]) { -// globTurbineData[iTurb].TurbID = turbNode["turb_id"].as(); -// } else { -// std::cerr << "turb_id not present for Turbine " << iTurb << " in the input file" << std::endl; -// } - -// if (turbNode["FAST_input_filename"]) { -// globTurbineData[iTurb].FASTInputFileName = turbNode["FAST_input_filename"].as(); -// } else { -// std::cerr << "FAST_input_filename not present for Turbine " << iTurb << " in the input file" << std::endl; -// } - -// if (turbNode["restart_filename"]) { -// globTurbineData[iTurb].FASTRestartFileName = turbNode["restart_filename"].as(); -// } else { -// std::cerr << " not present for Turbine " << iTurb << " in the input file" << std::endl; -// } - -// if (turbNode["turbine_base_pos"].IsSequence() ) { -// globTurbineData[iTurb].TurbineBasePos = turbNode["turbine_base_pos"].as >() ; -// } else { -// std::cerr << "turbine_base_pos not present for Turbine " << iTurb << " in the input file" << std::endl; -// } - -// if (turbNode["turbine_hub_pos"].IsSequence() ) { -// globTurbineData[iTurb].TurbineHubPos = turbNode["turbine_hub_pos"].as >() ; -// } else { -// std::cerr << "turbine_hub_pos not present for Turbine " << iTurb << " in the input file" << std::endl; -// } - -// if (turbNode["num_force_pts_blade"]) { -// globTurbineData[iTurb].numForcePtsBlade = turbNode["num_force_pts_blade"].as(); -// } else { -// std::cout << "num_force_pts_blade not present for Turbine " << iTurb << " in the input file" << std::endl; -// std::cout << " ... proceeding anyway. Will fail if you want CFD integration" << std::endl ; -// } - -// if (turbNode["num_force_pts_tower"]) { -// globTurbineData[iTurb].numForcePtsTwr = turbNode["num_force_pts_tower"].as(); -// } else { -// std::cout << "num_force_pts_tower not present for Turbine " << iTurb << " in the input file" << std::endl; -// std::cout << " ... proceeding anyway. Will fail if you want CFD integration" << std::endl ; -// } - -// return ; -// } - - void fast::OpenFAST::allocateTurbinesToProcsSimple() { // Allocate turbines to each processor - round robin fashion @@ -669,34 +621,34 @@ void fast::OpenFAST::allocateTurbinesToProcsSimple() { void fast::OpenFAST::end() { - // Deallocate types we allocated earlier + // Deallocate types we allocated earlier - if ( !dryRun) { - bool stopTheProgram = false; - for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { - FAST_End(&iTurb, &stopTheProgram); - } - } - - MPI_Group_free(&fastMPIGroup); - if (MPI_COMM_NULL != fastMPIComm) { - MPI_Comm_free(&fastMPIComm); + if ( !dryRun) { + bool stopTheProgram = false; + for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { + FAST_End(&iTurb, &stopTheProgram); } - MPI_Group_free(&worldMPIGroup); - - if(scStatus) { - - destroy_SuperController(sc) ; - - if(scLibHandle != NULL) { - // close the library - std::cout << "Closing library...\n"; - dlclose(scLibHandle); - } - + } + + MPI_Group_free(&fastMPIGroup); + if (MPI_COMM_NULL != fastMPIComm) { + MPI_Comm_free(&fastMPIComm); + } + MPI_Group_free(&worldMPIGroup); + + if(scStatus) { + + destroy_SuperController(sc) ; + + if(scLibHandle != NULL) { + // close the library + std::cout << "Closing library...\n"; + dlclose(scLibHandle); } - + } + +} void fast::OpenFAST::loadSuperController(const fast::fastInputs & fi) { From 8a70af88df9796f223d50b30a231110ae04f4e4d Mon Sep 17 00:00:00 2001 From: Ganesh Vijayakumar Date: Mon, 15 May 2017 11:35:47 -0600 Subject: [PATCH 3/8] Updating OpenFOAM module to get new units for file I/O using GetNewUnit() --- modules-local/openfoam/src/OpenFOAM.f90 | 40 ++++++++++++++++++------- 1 file changed, 29 insertions(+), 11 deletions(-) diff --git a/modules-local/openfoam/src/OpenFOAM.f90 b/modules-local/openfoam/src/OpenFOAM.f90 index a949bb571e..9f7cb9c80e 100644 --- a/modules-local/openfoam/src/OpenFOAM.f90 +++ b/modules-local/openfoam/src/OpenFOAM.f90 @@ -458,6 +458,9 @@ SUBROUTINE SetOpFMForces(p_FAST, p_AD14, u_AD14, y_AD14, u_AD, y_AD, y_ED, OpFM, INTEGER(IntKi) :: J ! Loops through nodes / elements INTEGER(IntKi) :: K ! Loops through blades. INTEGER(IntKi) :: Node ! Node number for blade/node on mesh +#ifdef DEBUG_OPENFOAM + INTEGER(IntKi) :: actForcesFile, aerodynForcesFile ! Unit numbers for files containing actuator forces and aerodyn forces +#endif INTEGER(IntKi) :: ErrStat2 ! temporary Error status of the operation CHARACTER(ErrMsgLen) :: ErrMsg2 ! temporary Error message if ErrStat /= ErrID_None @@ -476,17 +479,23 @@ SUBROUTINE SetOpFMForces(p_FAST, p_AD14, u_AD14, y_AD14, u_AD, y_AD, y_ED, OpFM, ! blade nodes !....................... - open(unit=987,file='fast_aerodyn_velocity_forces.csv') - write(987,*) '#x, y, z, fx, fy, fz' +#ifdef DEBUG_OPENFOAM + CALL GetNewUnit( aerodynForcesFile ) + open(unit=aerodynForcesFile,file='fast_aerodyn_velocity_forces.csv') + write(aerodynForcesFile,*) '#x, y, z, fx, fy, fz' - open(unit=988,file='fast_actuator_forces.csv') - write(988,*) '#x, y, z, fx, fy, fz' + CALL GetNewUnit( actForcesFile ) + open(unit=actForcesFile,file='fast_actuator_forces.csv') + write(actForcesFile,*) '#x, y, z, fx, fy, fz' +#endif DO K = 1,OpFM%p%NumBl +#ifdef DEBUG_OPENFOAM DO J = 1,u_AD%BladeMotion(k)%NNodes - write(987,*) u_AD%BladeMotion(k)%TranslationDisp(1,j) + u_AD%BladeMotion(k)%Position(1,j), ', ', u_AD%BladeMotion(k)%TranslationDisp(2,j) + u_AD%BladeMotion(k)%Position(2,j), ', ', u_AD%BladeMotion(k)%TranslationDisp(3,j) + u_AD%BladeMotion(k)%Position(3,j), ', ', OpFM%y%u(1 + (k-1)*u_AD%BladeMotion(k)%NNodes + j), ', ', OpFM%y%v(1 + (k-1)*u_AD%BladeMotion(k)%NNodes + j), ', ', OpFM%y%w(1 + (k-1)*u_AD%BladeMotion(k)%NNodes + j), ', ', y_AD%BladeLoad(k)%Force(1,j), ', ', y_AD%BladeLoad(k)%Force(2,j), ', ', y_AD%BladeLoad(k)%Force(2,j) + write(aerodynForcesFile,*) u_AD%BladeMotion(k)%TranslationDisp(1,j) + u_AD%BladeMotion(k)%Position(1,j), ', ', u_AD%BladeMotion(k)%TranslationDisp(2,j) + u_AD%BladeMotion(k)%Position(2,j), ', ', u_AD%BladeMotion(k)%TranslationDisp(3,j) + u_AD%BladeMotion(k)%Position(3,j), ', ', OpFM%y%u(1 + (k-1)*u_AD%BladeMotion(k)%NNodes + j), ', ', OpFM%y%v(1 + (k-1)*u_AD%BladeMotion(k)%NNodes + j), ', ', OpFM%y%w(1 + (k-1)*u_AD%BladeMotion(k)%NNodes + j), ', ', y_AD%BladeLoad(k)%Force(1,j), ', ', y_AD%BladeLoad(k)%Force(2,j), ', ', y_AD%BladeLoad(k)%Force(2,j) END DO +#endif call Transfer_Line2_to_Line2( y_AD%BladeLoad(k), OpFM%m%ActForceLoads(k), OpFM%m%Line2_to_Line2_Loads(k), ErrStat2, ErrMsg2, u_AD%BladeMotion(k), OpFM%m%ActForceMotions(k) ) call Transfer_Line2_to_Point( OpFM%m%ActForceLoads(k), OpFM%m%ActForceLoadsPoints(k), OpFM%m%Line2_to_Point_Loads(k), ErrStat2, ErrMsg2, OpFM%m%ActForceMotions(k), OpFM%m%ActForceMotionsPoints(k) ) @@ -499,8 +508,11 @@ SUBROUTINE SetOpFMForces(p_FAST, p_AD14, u_AD14, y_AD14, u_AD, y_AD, y_ED, OpFM, OpFM%u%momentx(Node) = OpFM%m%ActForceLoadsPoints(k)%Moment(1,j) OpFM%u%momenty(Node) = OpFM%m%ActForceLoadsPoints(k)%Moment(2,j) OpFM%u%momentz(Node) = OpFM%m%ActForceLoadsPoints(k)%Moment(3,j) - - write(988,*) OpFM%u%pxForce(Node), ', ', OpFM%u%pyForce(Node), ', ', OpFM%u%pzForce(Node), ', ', OpFM%u%fx(Node), ', ', OpFM%u%fy(Node), ', ', OpFM%u%fz(Node), ', ' + +#ifdef DEBUG_OPENFOAM + write(actForcesFile,*) OpFM%u%pxForce(Node), ', ', OpFM%u%pyForce(Node), ', ', OpFM%u%pzForce(Node), ', ', OpFM%u%fx(Node), ', ', OpFM%u%fy(Node), ', ', OpFM%u%fz(Node), ', ' +#endif + END DO END DO !K = 1,OpFM%p%NumBl @@ -512,9 +524,11 @@ SUBROUTINE SetOpFMForces(p_FAST, p_AD14, u_AD14, y_AD14, u_AD, y_AD, y_ED, OpFM, ! mesh mapping from line2 mesh to point mesh k = SIZE(u_AD%BladeMotion) + 1 +#ifdef DEBUG_OPENFOAM DO J = 1,u_AD%TowerMotion%NNodes - write(987,*) u_AD%TowerMotion%TranslationDisp(1,j) + u_AD%TowerMotion%Position(1,j), ', ', u_AD%TowerMotion%TranslationDisp(2,j) + u_AD%TowerMotion%Position(2,j), ', ', u_AD%TowerMotion%TranslationDisp(3,j) + u_AD%TowerMotion%Position(3,j), ', ', OpFM%y%u(1 + OpFM%p%NumBl*u_AD%BladeMotion(k)%NNodes + j), ', ', OpFM%y%v(1 + OpFM%p%NumBl*u_AD%BladeMotion(k)%NNodes + j), ', ', OpFM%y%w(1 + OpFM%p%NumBl*u_AD%BladeMotion(k)%NNodes + j), ', ', y_AD%TowerLoad%Force(1,j), ', ', y_AD%TowerLoad%Force(2,j), ', ', y_AD%TowerLoad%Force(2,j) + write(aerodynForcesFile,*) u_AD%TowerMotion%TranslationDisp(1,j) + u_AD%TowerMotion%Position(1,j), ', ', u_AD%TowerMotion%TranslationDisp(2,j) + u_AD%TowerMotion%Position(2,j), ', ', u_AD%TowerMotion%TranslationDisp(3,j) + u_AD%TowerMotion%Position(3,j), ', ', OpFM%y%u(1 + OpFM%p%NumBl*u_AD%BladeMotion(k)%NNodes + j), ', ', OpFM%y%v(1 + OpFM%p%NumBl*u_AD%BladeMotion(k)%NNodes + j), ', ', OpFM%y%w(1 + OpFM%p%NumBl*u_AD%BladeMotion(k)%NNodes + j), ', ', y_AD%TowerLoad%Force(1,j), ', ', y_AD%TowerLoad%Force(2,j), ', ', y_AD%TowerLoad%Force(2,j) END DO +#endif call Transfer_Line2_to_Line2( y_AD%TowerLoad, OpFM%m%ActForceLoads(k), OpFM%m%Line2_to_Line2_Loads(k), ErrStat2, ErrMsg2, u_AD%TowerMotion, OpFM%m%ActForceMotions(k) ) call Transfer_Line2_to_Point( OpFM%m%ActForceLoads(k), OpFM%m%ActForceLoadsPoints(k), OpFM%m%Line2_to_Point_Loads(k), ErrStat2, ErrMsg2, OpFM%m%ActForceMotions(k), OpFM%m%ActForceMotionsPoints(k) ) @@ -528,11 +542,15 @@ SUBROUTINE SetOpFMForces(p_FAST, p_AD14, u_AD14, y_AD14, u_AD, y_AD, y_ED, OpFM, OpFM%u%momenty(Node) = OpFM%m%ActForceLoadsPoints(k)%Moment(2,j) OpFM%u%momentz(Node) = OpFM%m%ActForceLoadsPoints(k)%Moment(3,j) - write(988,*) OpFM%u%pxForce(Node), ', ', OpFM%u%pyForce(Node), ', ', OpFM%u%pzForce(Node), ', ', OpFM%u%fx(Node), ', ', OpFM%u%fy(Node), ', ', OpFM%u%fz(Node), ', ' +#ifdef DEBUG_OPENFOAM + write(actForcesFile,*) OpFM%u%pxForce(Node), ', ', OpFM%u%pyForce(Node), ', ', OpFM%u%pzForce(Node), ', ', OpFM%u%fx(Node), ', ', OpFM%u%fy(Node), ', ', OpFM%u%fz(Node), ', ' +#endif END DO - close(987) - close(988) +#ifdef DEBUG_OPENFOAM + close(aerodynForcesFile) + close(actForcesFile) +#endif END SUBROUTINE SetOpFMForces !---------------------------------------------------------------------------------------------------------------------------------- From 9082e5c63c3705d081d33d1bec061f33cc93e3bd Mon Sep 17 00:00:00 2001 From: Ganesh Vijayakumar Date: Tue, 16 May 2017 17:53:40 -0600 Subject: [PATCH 4/8] Updating the debug write out of processors containing non-zero instances of FAST --- glue-codes/fast-cpp/src/OpenFAST.cpp | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/glue-codes/fast-cpp/src/OpenFAST.cpp b/glue-codes/fast-cpp/src/OpenFAST.cpp index 602d3ed569..88cc218ee9 100644 --- a/glue-codes/fast-cpp/src/OpenFAST.cpp +++ b/glue-codes/fast-cpp/src/OpenFAST.cpp @@ -546,23 +546,25 @@ void fast::OpenFAST::allocateMemory() { int nProcsWithTurbines=0; turbineProcs.resize(turbineSetProcs.size()); - std::ofstream turbineAllocFile; - turbineAllocFile.open("turbineAlloc." + std::to_string(worldMPIRank) + ".txt") ; + for (std::set::const_iterator p = turbineSetProcs.begin(); p != turbineSetProcs.end(); p++) { turbineProcs[nProcsWithTurbines] = *p; + nProcsWithTurbines++ ; + } - if (dryRun) { - if ( worldMPIRank == turbineProcs[nProcsWithTurbines] ) { - for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { - turbineAllocFile << "Proc " << worldMPIRank << " loc iTurb " << iTurb << " glob iTurb " << turbineMapProcToGlob[iTurb] << std::endl ; - } + if (dryRun) { + if (nTurbinesProc > 0) { + std::ofstream turbineAllocFile; + turbineAllocFile.open("turbineAlloc." + std::to_string(worldMPIRank) + ".txt") ; + for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { + turbineAllocFile << "Proc " << worldMPIRank << " loc iTurb " << iTurb << " glob iTurb " << turbineMapProcToGlob[iTurb] << std::endl ; } + turbineAllocFile.flush(); + turbineAllocFile.close() ; } - - nProcsWithTurbines++ ; + } - turbineAllocFile.flush(); - turbineAllocFile.close() ; + // Construct a group containing all procs running atleast 1 turbine in FAST MPI_Group_incl(worldMPIGroup, nProcsWithTurbines, &turbineProcs[0], &fastMPIGroup) ; From 9094d734b3250ecfedc0d86cbe508b3a8065bd2b Mon Sep 17 00:00:00 2001 From: Ganesh Vijayakumar Date: Fri, 19 May 2017 14:17:31 -0600 Subject: [PATCH 5/8] Adding an option to set velocity directly at the velocity nodes or at the force nodes --- glue-codes/fast-cpp/src/OpenFAST.H | 1 + glue-codes/fast-cpp/src/OpenFAST.cpp | 10 ++++++++++ 2 files changed, 11 insertions(+) diff --git a/glue-codes/fast-cpp/src/OpenFAST.H b/glue-codes/fast-cpp/src/OpenFAST.H index 709fb0bee7..5107c68c82 100644 --- a/glue-codes/fast-cpp/src/OpenFAST.H +++ b/glue-codes/fast-cpp/src/OpenFAST.H @@ -158,6 +158,7 @@ class OpenFAST { ActuatorNodeType getVelNodeType(int iTurbGlob, int iNode); void getVelNodeCoordinates(std::vector & currentCoords, int iNode, int iTurbGlob); void setVelocity(std::vector & velocity, int iNode, int iTurbGlob); + void setVelocityForceNode(std::vector & velocity, int iNode, int iTurbGlob); void interpolateVel_ForceToVelNodes(); ActuatorNodeType getForceNodeType(int iTurbGlob, int iNode); void getForceNodeCoordinates(std::vector & currentCoords, int iNode, int iTurbGlob); diff --git a/glue-codes/fast-cpp/src/OpenFAST.cpp b/glue-codes/fast-cpp/src/OpenFAST.cpp index 88cc218ee9..d728cd05bd 100644 --- a/glue-codes/fast-cpp/src/OpenFAST.cpp +++ b/glue-codes/fast-cpp/src/OpenFAST.cpp @@ -339,6 +339,16 @@ void fast::OpenFAST::setVelocity(std::vector & currentVelocity, int iNod // Set velocity at current node of current turbine - int iTurbLoc = get_localTurbNo(iTurbGlob); for(int j=0; j < iTurbLoc; j++) iNode = iNode - get_numVelPtsLoc(iTurbLoc); + cDriver_Output_to_FAST[iTurbLoc].u[iNode] = currentVelocity[0]; + cDriver_Output_to_FAST[iTurbLoc].v[iNode] = currentVelocity[1]; + cDriver_Output_to_FAST[iTurbLoc].w[iNode] = currentVelocity[2]; +} + +void fast::OpenFAST::setVelocityForceNode(std::vector & currentVelocity, int iNode, int iTurbGlob) { + + // Set velocity at current node of current turbine - + int iTurbLoc = get_localTurbNo(iTurbGlob); + for(int j=0; j < iTurbLoc; j++) iNode = iNode - get_numForcePtsLoc(iTurbLoc); forceNodeVel[iTurbLoc][iNode][0] = currentVelocity[0]; forceNodeVel[iTurbLoc][iNode][1] = currentVelocity[1]; forceNodeVel[iTurbLoc][iNode][2] = currentVelocity[2]; From 135c751cd2c2d2f65f46543eca817f1c3e3bc35d Mon Sep 17 00:00:00 2001 From: Ganesh Vijayakumar Date: Mon, 22 May 2017 07:23:45 -0600 Subject: [PATCH 6/8] Reverting to old style std vector declaration and adding some ifndefs to allow for compilation with older compilers --- glue-codes/fast-cpp/src/OpenFAST.H | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/glue-codes/fast-cpp/src/OpenFAST.H b/glue-codes/fast-cpp/src/OpenFAST.H index 5107c68c82..04504dcc31 100644 --- a/glue-codes/fast-cpp/src/OpenFAST.H +++ b/glue-codes/fast-cpp/src/OpenFAST.H @@ -12,8 +12,12 @@ #include #include #include "dlfcn.h" -#define OMPI_SKIP_MPICXX -#define MPICH_SKIP_MPICXX +#ifndef OMPI_SKIP_MPICXX + #define OMPI_SKIP_MPICXX +#endif +#ifndef MPICH_SKIP_MPICXX + #define MPICH_SKIP_MPICXX +#endif #include "mpi.h" #include "SC.h" @@ -79,8 +83,8 @@ class OpenFAST { bool timeZero; double dtFAST; double tMax; - std::vector> TurbineBasePos; - std::vector> TurbineHubPos; + std::vector > TurbineBasePos; + std::vector > TurbineHubPos; std::vector TurbID; std::vector FASTInputFileName; std::vector CheckpointFileRoot; @@ -98,7 +102,7 @@ class OpenFAST { std::vector scOutputsGlob; // # outputs from the supercontroller for all turbines std::vector scInputsGlob; // # inputs to the supercontroller for all turbines - std::vector>> forceNodeVel; // Velocity at force nodes - Store temporarily to interpolate to the velocity nodes + std::vector > > forceNodeVel; // Velocity at force nodes - Store temporarily to interpolate to the velocity nodes std::vector cDriver_Input_from_FAST; std::vector cDriver_Output_to_FAST; From b9b5e626e5fe855b03d9a0be901dad6b850cfe16 Mon Sep 17 00:00:00 2001 From: Ganesh Vijayakumar Date: Mon, 22 May 2017 18:12:58 -0600 Subject: [PATCH 7/8] Another bugfix in sending orientations at force nodes back to the driver code --- glue-codes/fast-cpp/src/OpenFAST.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/glue-codes/fast-cpp/src/OpenFAST.cpp b/glue-codes/fast-cpp/src/OpenFAST.cpp index d728cd05bd..9911005c9f 100644 --- a/glue-codes/fast-cpp/src/OpenFAST.cpp +++ b/glue-codes/fast-cpp/src/OpenFAST.cpp @@ -309,7 +309,7 @@ void fast::OpenFAST::getForceNodeOrientation(std::vector & currentOrient int iTurbLoc = get_localTurbNo(iTurbGlob); for(int j=0; j < iTurbLoc; j++) iNode = iNode - get_numForcePtsLoc(iTurbLoc); for(int i=0;i<9;i++) { - currentOrientation[i] = cDriver_Input_from_FAST[iTurbLoc].pxForce[iNode*9+i] ; + currentOrientation[i] = cDriver_Input_from_FAST[iTurbLoc].pOrientation[iNode*9+i] ; } } From 15d9d20349b265aac11b2b9fb3cd36e3129e6722 Mon Sep 17 00:00:00 2001 From: Ganesh Vijayakumar Date: Tue, 23 May 2017 10:57:15 -0600 Subject: [PATCH 8/8] Allowing for new restart capability in OpenFAST - C++ API While OpenFAST does have a restart capability, some of our work is dependent on the use of a third party turbine controller using the Bladed interface. The Bladed interface could have a restart capability, but is not guaranteed to have one. To work around this limitation, the OpenFAST-C++ API will provide a new capability as shown below: When the C++ driver starts from time zero, OpenFAST will start normally and write the velocity data the Aerodyn nodes to a file every time step. When the C++ driver restarts from a non-zero time (time > 0), we have two options 1. trueRestart - All components of FAST including the turbine controller are expected to have a checkpoint file at the specified start time. The C++ driver will call OpenFAST to restart from the checkpoint files at the start time and then proceed to step together with OpenFAST for the rest of the simulation. 2. restartDriverInitFAST - The C++ driver will call OpenFAST to start from scratch (time = 0), read the file containing the velocity data at the Aerodyn nodes, and then step OpenFAST upto the start time of the simulation by applying the velocity data at the Aerodyn nodes. Once OpenFAST reaches the start time of the simulation, the C++ driver and OpenFAST will continue to step together normally while continuing to write out the velocity data the Aerodyn nodes at every time step. --- glue-codes/fast-cpp/CMakeLists.txt | 5 +- glue-codes/fast-cpp/src/FAST_Prog.cpp | 14 +- glue-codes/fast-cpp/src/OpenFAST.H | 23 ++- glue-codes/fast-cpp/src/OpenFAST.cpp | 277 ++++++++++++++++++++++++-- 4 files changed, 297 insertions(+), 22 deletions(-) diff --git a/glue-codes/fast-cpp/CMakeLists.txt b/glue-codes/fast-cpp/CMakeLists.txt index 613a0ecf00..8d75a85877 100644 --- a/glue-codes/fast-cpp/CMakeLists.txt +++ b/glue-codes/fast-cpp/CMakeLists.txt @@ -18,7 +18,7 @@ set(CMAKE_CXX_STANDARD 11) set(CMAKE_CXX_STANDARD_REQUIRED ON) find_package(MPI REQUIRED) -find_package(HDF5 REQUIRED) +find_package(HDF5 REQUIRED PATHS HDF5_ROOT) find_package(YAMLCPP REQUIRED) include_directories(${YAML_INCLUDES}) @@ -32,7 +32,8 @@ add_library(openfastcpplib src/OpenFAST.cpp) target_link_libraries(openfastcpplib openfastlib - ${HDF5_LIBRARIES} + ${HDF5_C_STATIC_LIBRARY} + ${HDF5_HL_STATIC_LIBRARY} ${CMAKE_DL_LIBS}) add_executable(openfastcpp diff --git a/glue-codes/fast-cpp/src/FAST_Prog.cpp b/glue-codes/fast-cpp/src/FAST_Prog.cpp index 06993d8509..2b6b79394d 100644 --- a/glue-codes/fast-cpp/src/FAST_Prog.cpp +++ b/glue-codes/fast-cpp/src/FAST_Prog.cpp @@ -45,6 +45,18 @@ void readInputFile(fast::fastInputs & fi, std::string cInterfaceInputFile, doubl if(cDriverInp["debug"]) { fi.debug = cDriverInp["debug"].as(); } + + if(cDriverInp["simStart"]) { + if (cDriverInp["simStart"].as() == "init") { + fi.simStart = fast::init; + } else if(cDriverInp["simStart"].as() == "trueRestart") { + fi.simStart = fast::trueRestart; + } else if(cDriverInp["simStart"].as() == "restartDriverInitFAST") { + fi.simStart = fast::restartDriverInitFAST; + } else { + throw std::runtime_error("simStart is not well defined in the input file"); + } + } fi.tStart = cDriverInp["tStart"].as(); *tEnd = cDriverInp["tEnd"].as(); @@ -110,7 +122,7 @@ int main() { // Or allocate turbines to procs by calling "setTurbineProcNo(iTurbGlob, procId)" for turbine. FAST.init(); - if (!FAST.isRestart()) { + if (FAST.isTimeZero()) { FAST.solution0(); } diff --git a/glue-codes/fast-cpp/src/OpenFAST.H b/glue-codes/fast-cpp/src/OpenFAST.H index 04504dcc31..e9a605cc43 100644 --- a/glue-codes/fast-cpp/src/OpenFAST.H +++ b/glue-codes/fast-cpp/src/OpenFAST.H @@ -41,6 +41,14 @@ enum ActuatorNodeType { ActuatorNodeType_END }; +enum simStartType { + init = 0, + trueRestart = 1, + restartDriverInitFAST = 2, + simStartType_END +}; + + class fastInputs { public: @@ -50,6 +58,7 @@ class fastInputs { bool dryRun; bool debug; double tStart; + simStartType simStart; int nEveryCheckPoint; double tMax; double dtFAST; @@ -79,7 +88,7 @@ class OpenFAST { std::vector globTurbineData; int nTurbinesProc; int nTurbinesGlob; - bool restart; + simStartType simStart; bool timeZero; double dtFAST; double tMax; @@ -103,6 +112,9 @@ class OpenFAST { std::vector scInputsGlob; // # inputs to the supercontroller for all turbines std::vector > > forceNodeVel; // Velocity at force nodes - Store temporarily to interpolate to the velocity nodes + std::vector > velNodeData; // Position and velocity data at the velocity (aerodyn) nodes - (nTurbines, nTimesteps * nPoints * 6) + hid_t velNodeDataFile; // HDF-5 tag of file containing velocity (aerodyn) node data file + std::vector cDriver_Input_from_FAST; std::vector cDriver_Output_to_FAST; @@ -151,8 +163,14 @@ class OpenFAST { void init(); void solution0(); void step(); + void stepNoWrite(); void end(); + hid_t openVelocityDataFile(bool createFile); + int readVelocityData(); + void writeVelocityData(hid_t h5file, int iTurb, int iTimestep, OpFM_InputType_t iData, OpFM_OutputType_t oData); + herr_t closeVelocityDataFile(int nt_global, hid_t velDataFile); + void setTurbineProcNo(int iTurbGlob, int procNo) { turbineMapGlobToProc[iTurbGlob] = procNo; } void allocateTurbinesToProcsSimple(); void getApproxHubPos(std::vector & currentCoords, int iTurbGlob); @@ -173,7 +191,7 @@ class OpenFAST { int get_ntStart() { return ntStart; } bool isDryRun() { return dryRun; } bool isDebug() { return debug; } - bool isRestart() { return restart; } + simStartType get_simStartType() { return simStart; } bool isTimeZero() { return timeZero; } int get_procNo(int iTurbGlob) { return turbineMapGlobToProc[iTurbGlob] ; } // Get processor number of a turbine with global id 'iTurbGlob' int get_localTurbNo(int iTurbGlob) { return reverseTurbineMapProcToGlob[iTurbGlob]; } @@ -209,6 +227,7 @@ class OpenFAST { void fillScOutputsLoc() ; void setOutputsToFAST(OpFM_InputType_t cDriver_Input_from_FAST, OpFM_OutputType_t cDriver_Output_to_FAST) ; // An example to set velocities at the Aerodyn nodes + void applyVelocityData(int iPrestart, int iTurb, OpFM_OutputType_t cDriver_Output_to_FAST, std::vector & velData) ; }; diff --git a/glue-codes/fast-cpp/src/OpenFAST.cpp b/glue-codes/fast-cpp/src/OpenFAST.cpp index 9911005c9f..2e18fa845a 100644 --- a/glue-codes/fast-cpp/src/OpenFAST.cpp +++ b/glue-codes/fast-cpp/src/OpenFAST.cpp @@ -23,7 +23,7 @@ fast::OpenFAST::OpenFAST(): nTurbinesGlob(0), nTurbinesProc(0), scStatus(false), -restart(false), +simStart(fast::init), timeZero(false) { } @@ -38,8 +38,9 @@ void fast::OpenFAST::init() { allocateMemory(); if (!dryRun) { - // If restart - if (restart == true) { + switch (simStart) { + + case fast::trueRestart: for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { /* note that this will set nt_global inside the FAST library */ @@ -50,17 +51,46 @@ void fast::OpenFAST::init() { int nfpts = get_numForcePtsLoc(iTurb); forceNodeVel[iTurb].resize(nfpts); for (int k = 0; k < nfpts; k++) forceNodeVel[iTurb][k].resize(3) ; + } + if (nTurbinesProc > 0) velNodeDataFile = openVelocityDataFile(false); + if(scStatus) { sc->readRestartFile(nt_global); } - + + break ; - } else { + case fast::init: // this calls the Init() routines of each module + for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { + FAST_OpFM_Init(&iTurb, &tMax, FASTInputFileName[iTurb].data(), &TurbID[iTurb], &numScOutputs, &numScInputs, &numForcePtsBlade[iTurb], &numForcePtsTwr[iTurb], TurbineBasePos[iTurb].data(), &AbortErrLev, &dtFAST, &numBlades[iTurb], &numVelPtsBlade[iTurb], &cDriver_Input_from_FAST[iTurb], &cDriver_Output_to_FAST[iTurb], &cDriverSC_Input_from_FAST[iTurb], &cDriverSC_Output_to_FAST[iTurb], &ErrStat, ErrMsg); + checkError(ErrStat, ErrMsg); + + timeZero = true; + + numVelPtsTwr[iTurb] = cDriver_Output_to_FAST[iTurb].u_Len - numBlades[iTurb]*numVelPtsBlade[iTurb] - 1; + + int nfpts = get_numForcePtsLoc(iTurb); + forceNodeVel[iTurb].resize(nfpts); + for (int k = 0; k < nfpts; k++) forceNodeVel[iTurb][k].resize(3) ; + + if ( isDebug() ) { + for (int iNode=0; iNode < get_numVelPtsLoc(iTurb); iNode++) { + std::cout << "Node " << iNode << " Position = " << cDriver_Input_from_FAST[iTurb].pxVel[iNode] << " " << cDriver_Input_from_FAST[iTurb].pyVel[iNode] << " " << cDriver_Input_from_FAST[iTurb].pzVel[iNode] << " " << std::endl ; + } + } + } + + if (nTurbinesProc > 0) velNodeDataFile = openVelocityDataFile(true); + + break ; + + case fast::restartDriverInitFAST: + for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { FAST_OpFM_Init(&iTurb, &tMax, FASTInputFileName[iTurb].data(), &TurbID[iTurb], &numScOutputs, &numScInputs, &numForcePtsBlade[iTurb], &numForcePtsTwr[iTurb], TurbineBasePos[iTurb].data(), &AbortErrLev, &dtFAST, &numBlades[iTurb], &numVelPtsBlade[iTurb], &cDriver_Input_from_FAST[iTurb], &cDriver_Output_to_FAST[iTurb], &cDriverSC_Input_from_FAST[iTurb], &cDriverSC_Output_to_FAST[iTurb], &ErrStat, ErrMsg); checkError(ErrStat, ErrMsg); @@ -79,10 +109,37 @@ void fast::OpenFAST::init() { } } } + + int nTimesteps; - } - } + if (nTurbinesProc > 0) { + nTimesteps = readVelocityData(); + if (nTimesteps != ntStart) throw std::runtime_error("Start time is not consistent with the number of timesteps in the velocity data file"); + } + for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { + applyVelocityData(0, iTurb, cDriver_Output_to_FAST[iTurb], velNodeData[iTurb]); + } + solution0() ; + + for (int iPrestart=0 ; iPrestart < ntStart; iPrestart++) { + for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { + applyVelocityData(iPrestart, iTurb, cDriver_Output_to_FAST[iTurb], velNodeData[iTurb]); + } + stepNoWrite(); + } + + if (nTurbinesProc > 0) velNodeDataFile = openVelocityDataFile(false); + + break; + + case fast::simStartType_END: + + break; + + } + + } } void fast::OpenFAST::solution0() { @@ -132,9 +189,12 @@ void fast::OpenFAST::step() { // set wind speeds at original locations // setOutputsToFAST(cDriver_Input_from_FAST[iTurb], cDriver_Output_to_FAST[iTurb]); - + // this advances the states, calls CalcOutput, and solves for next inputs. Predictor-corrector loop is imbeded here: // (note OpenFOAM could do subcycling around this step) + + writeVelocityData(velNodeDataFile, iTurb, nt_global, cDriver_Input_from_FAST[iTurb], cDriver_Output_to_FAST[iTurb]); + if ( isDebug() ) { std::ofstream fastcpp_velocity_file; @@ -185,6 +245,38 @@ void fast::OpenFAST::step() { } +void fast::OpenFAST::stepNoWrite() { + + /* ****************************** + set inputs from this code and call FAST: + ********************************* */ + + if(scStatus) { + sc->calcOutputs(scOutputsGlob); + fillScOutputsLoc(); + } + + for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { + + // set wind speeds at original locations + // setOutputsToFAST(cDriver_Input_from_FAST[iTurb], cDriver_Output_to_FAST[iTurb]); + + // this advances the states, calls CalcOutput, and solves for next inputs. Predictor-corrector loop is imbeded here: + // (note OpenFOAM could do subcycling around this step) + FAST_OpFM_Step(&iTurb, &ErrStat, ErrMsg); + checkError(ErrStat, ErrMsg); + + } + + if(scStatus) { + sc->updateStates(scInputsGlob); // Go from 'n' to 'n+1' based on input at previous time step + fillScInputsGlob(); // Update inputs to super controller for 'n+1' + } + + nt_global = nt_global + 1; + +} + void fast::OpenFAST::setInputs(const fast::fastInputs & fi ) { @@ -202,6 +294,7 @@ void fast::OpenFAST::setInputs(const fast::fastInputs & fi ) { debug = fi.debug; tStart = fi.tStart; + simStart = fi.simStart; nEveryCheckPoint = fi.nEveryCheckPoint; tMax = fi.tMax; loadSuperController(fi); @@ -209,11 +302,12 @@ void fast::OpenFAST::setInputs(const fast::fastInputs & fi ) { ntStart = int(tStart/dtFAST); - nt_global = ntStart; - if (ntStart > 0) { - restart = true ; - } - + if (simStart == fast::restartDriverInitFAST) { + nt_global = 0; + } else { + nt_global = ntStart; + } + globTurbineData.resize(nTurbinesGlob); globTurbineData = fi.globTurbineData; @@ -245,10 +339,10 @@ void fast::OpenFAST::setOutputsToFAST(OpFM_InputType_t cDriver_Input_from_FAST, cDriver_Output_to_FAST.w[j] = 0.0; } - // call supercontroller - for (int j = 0; j < cDriver_Output_to_FAST.SuperController_Len; j++){ - cDriver_Output_to_FAST.SuperController[j] = (float) j; // set it somehow.... (would be set from the SuperController outputs) - } + // // call supercontroller + // for (int j = 0; j < cDriver_Output_to_FAST.SuperController_Len; j++){ + // cDriver_Output_to_FAST.SuperController[j] = (float) j; // set it somehow.... (would be set from the SuperController outputs) + // } } @@ -635,6 +729,8 @@ void fast::OpenFAST::end() { // Deallocate types we allocated earlier + if (nTurbinesProc > 0) closeVelocityDataFile(nt_global, velNodeDataFile); + if ( !dryRun) { bool stopTheProgram = false; for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { @@ -662,6 +758,153 @@ void fast::OpenFAST::end() { } +int fast::OpenFAST::readVelocityData() { + + int nTurbines, nTimesteps; + + hid_t velDataFile = H5Fopen(("velDatafile." + std::to_string(worldMPIRank) + ".h5").c_str(), H5F_ACC_RDWR, H5P_DEFAULT); + + { + hid_t attr = H5Aopen(velDataFile, "nTurbines", H5P_DEFAULT); + herr_t ret = H5Aread(attr, H5T_NATIVE_INT, &nTurbines) ; + H5Aclose(attr); + + attr = H5Aopen(velDataFile, "nTimesteps", H5P_DEFAULT); + ret = H5Aread(attr, H5T_NATIVE_INT, &nTimesteps) ; + H5Aclose(attr); + + } + + // Allocate memory and read the velocity data. + velNodeData.resize(nTurbines); + for (int iTurb=0; iTurb < nTurbines; iTurb++) { + int nVelPts = get_numVelPtsLoc(iTurb) ; + velNodeData[iTurb].resize(nTimesteps*nVelPts*6) ; + hid_t dset_id = H5Dopen2(velDataFile, ("/turbine" + std::to_string(iTurb)).c_str(), H5P_DEFAULT); + hid_t dspace_id = H5Dget_space(dset_id); + + hsize_t start[3]; start[1] = 0; start[2] = 0; + hsize_t count[3]; count[0] = 1; count[1] = nVelPts; count[2] = 6; + hid_t mspace_id = H5Screate_simple(3, count, NULL); + + for (int iStep=0; iStep < nTimesteps; iStep++) { + start[0] = iStep; + H5Sselect_hyperslab(dspace_id, H5S_SELECT_SET, start, NULL, count, NULL); + herr_t status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, mspace_id, dspace_id, H5P_DEFAULT, &velNodeData[iTurb][iStep*nVelPts*6] ); + } + herr_t status = H5Dclose(dset_id); + + + } + + return nTimesteps; + +} + +hid_t fast::OpenFAST::openVelocityDataFile(bool createFile) { + + hid_t velDataFile; + if (createFile) { + // Open the file in create mode + velDataFile = H5Fcreate(("velDatafile." + std::to_string(worldMPIRank) + ".h5").c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + + { + hsize_t dims[1]; + dims[0] = 1; + hid_t dataSpace = H5Screate_simple(1, dims, NULL); + hid_t attr = H5Acreate2(velDataFile, "nTurbines", H5T_NATIVE_INT, dataSpace, H5P_DEFAULT, H5P_DEFAULT) ; + herr_t status = H5Awrite(attr, H5T_NATIVE_INT, &nTurbinesProc); + status = H5Aclose(attr); + status = H5Sclose(dataSpace); + + dataSpace = H5Screate_simple(1, dims, NULL); + attr = H5Acreate2(velDataFile, "nTimesteps", H5T_NATIVE_INT, dataSpace, H5P_DEFAULT, H5P_DEFAULT) ; + status = H5Aclose(attr); + status = H5Sclose(dataSpace); + } + + int ntMax = tMax/dtFAST ; + + for (int iTurb = 0; iTurb < nTurbinesProc; iTurb++) { + int nVelPts = get_numVelPtsLoc(iTurb); + hsize_t dims[3]; + dims[0] = ntMax; dims[1] = nVelPts; dims[2] = 6 ; + + hsize_t chunk_dims[3]; + chunk_dims[0] = 1; chunk_dims[1] = nVelPts; chunk_dims[2] = 6; + hid_t dcpl_id = H5Pcreate(H5P_DATASET_CREATE); + H5Pset_chunk(dcpl_id, 3, chunk_dims); + + hid_t dataSpace = H5Screate_simple(3, dims, NULL); + hid_t dataSet = H5Dcreate(velDataFile, ("/turbine" + std::to_string(iTurb)).c_str(), H5T_NATIVE_DOUBLE, dataSpace, H5P_DEFAULT, dcpl_id, H5P_DEFAULT); + + herr_t status = H5Pclose(dcpl_id); + status = H5Dclose(dataSet); + status = H5Sclose(dataSpace); + } + + } else { + // Open the file in append mode + velDataFile = H5Fopen(("velDatafile." + std::to_string(worldMPIRank) + ".h5").c_str(), H5F_ACC_RDWR, H5P_DEFAULT); + } + + return velDataFile; + +} + +herr_t fast::OpenFAST::closeVelocityDataFile(int nt_global, hid_t velDataFile) { + + + hid_t attr_id = H5Aopen_by_name(velDataFile, ".", "nTimesteps", H5P_DEFAULT, H5P_DEFAULT); + herr_t status = H5Awrite(attr_id, H5T_NATIVE_INT, &nt_global); + status = H5Aclose(attr_id); + + status = H5Fclose(velDataFile) ; + return status; +} + + +void fast::OpenFAST::writeVelocityData(hid_t h5File, int iTurb, int iTimestep, OpFM_InputType_t iData, OpFM_OutputType_t oData) { + + + hsize_t start[3]; start[0] = iTimestep; start[1] = 0; start[2] = 0; + int nVelPts = get_numVelPtsLoc(iTurb) ; + hsize_t count[3]; count[0] = 1; count[1] = nVelPts; count[2] = 6; + + std::vector tmpVelData; + tmpVelData.resize(nVelPts * 6); + + for (int iNode=0 ; iNode < nVelPts; iNode++) { + tmpVelData[iNode*6 + 0] = iData.pxVel[iNode]; + tmpVelData[iNode*6 + 1] = iData.pyVel[iNode]; + tmpVelData[iNode*6 + 2] = iData.pzVel[iNode]; + tmpVelData[iNode*6 + 3] = oData.u[iNode]; + tmpVelData[iNode*6 + 4] = oData.v[iNode]; + tmpVelData[iNode*6 + 5] = oData.w[iNode]; + } + + hid_t dset_id = H5Dopen2(h5File, ("/turbine" + std::to_string(iTurb)).c_str(), H5P_DEFAULT); + hid_t dspace_id = H5Dget_space(dset_id); + H5Sselect_hyperslab(dspace_id, H5S_SELECT_SET, start, NULL, count, NULL); + hid_t mspace_id = H5Screate_simple(3, count, NULL); + H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, mspace_id, dspace_id, H5P_DEFAULT, tmpVelData.data()); + + H5Dclose(dset_id); + H5Sclose(dspace_id); + H5Sclose(mspace_id); + +} + +void fast::OpenFAST::applyVelocityData(int iPrestart, int iTurb, OpFM_OutputType_t cDriver_Output_to_FAST, std::vector & velData) { + + int nVelPts = get_numVelPtsLoc(iTurb); + for (int j = 0; j < nVelPts; j++){ + cDriver_Output_to_FAST.u[j] = velData[(iPrestart*nVelPts+j)*6 + 3]; + cDriver_Output_to_FAST.v[j] = velData[(iPrestart*nVelPts+j)*6 + 4]; + cDriver_Output_to_FAST.w[j] = velData[(iPrestart*nVelPts+j)*6 + 5]; + } + +} void fast::OpenFAST::loadSuperController(const fast::fastInputs & fi) {