diff --git a/.github/workflows/automated-dev-tests.yml b/.github/workflows/automated-dev-tests.yml index 1a7c6f7baa..7d8b46052c 100644 --- a/.github/workflows/automated-dev-tests.yml +++ b/.github/workflows/automated-dev-tests.yml @@ -146,3 +146,56 @@ jobs: - name: Test working-directory: ${{runner.workspace}}/build run: ./glue-codes/openfast/openfast -v + + cpp-interface-tests: + runs-on: ubuntu-20.04 + steps: + - name: Checkout + uses: actions/checkout@main + with: + submodules: recursive + + - name: Setup Python + uses: actions/setup-python@v2 + with: + python-version: '3.7' + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install numpy Bokeh==1.4 + sudo apt-get install libhdf5-dev libopenmpi-dev libyaml-cpp-dev + + - name: Setup Workspace + run: cmake -E make_directory ${{runner.workspace}}/build + - name: Configure Build + working-directory: ${{runner.workspace}}/build + run: | + cmake \ + -DCMAKE_INSTALL_PREFIX:PATH=${{runner.workspace}}/install \ + -DCMAKE_Fortran_COMPILER:STRING=${{env.FORTRAN_COMPILER}} \ + -DCMAKE_BUILD_TYPE:STRING=RelWithDebInfo \ + -DBUILD_OPENFAST_CPP_API:BOOL=ON \ + -DBUILD_SHARED_LIBS:BOOL=ON \ + -DBUILD_TESTING:BOOL=ON \ + -DCTEST_PLOT_ERRORS:BOOL=ON \ + ${GITHUB_WORKSPACE} + - name: Build OpenFAST C++ API + # if: contains(github.event.head_commit.message, 'Action - Test All') || contains(github.event.pull_request.labels.*.name, 'Action - Test All') + working-directory: ${{runner.workspace}}/build + run: | + cmake --build . --target openfastcpp -- -j ${{env.NUM_PROCS}} + cmake --build . --target regression_tests -- -j ${{env.NUM_PROCS}} + + - name: Run OpenFAST C++ API tests + working-directory: ${{runner.workspace}}/build + run: | + ctest -VV -L cpp + + - name: Failing test artifacts + uses: actions/upload-artifact@v2 + if: failure() + with: + name: test-results + path: | + ${{runner.workspace}}/build/reg_tests/glue-codes/openfast-cpp + !${{runner.workspace}}/build/reg_tests/glue-codes/openfast-cpp/5MW_Baseline diff --git a/cmake/FindYAMLCPP.cmake b/cmake/FindYAMLCPP.cmake deleted file mode 100644 index a3557319e3..0000000000 --- a/cmake/FindYAMLCPP.cmake +++ /dev/null @@ -1,27 +0,0 @@ -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# - -find_path(YAML_INCLUDES - yaml-cpp/yaml.h - HINTS ${YAML_ROOT} ${CMAKE_INSTALL_PREFIX} - PATH_SUFFIXES include) - -find_library(YAML_LIBRARIES - NAMES yaml-cpp libyaml-cpp.a - HINTS ${YAML_ROOT} ${CMAKE_INSTALL_PREFIX} - PATH_SUFFIXES lib) - -include(FindPackageHandleStandardArgs) -find_package_handle_standard_args(YAMLCPP DEFAULT_MSG YAML_INCLUDES YAML_LIBRARIES) -mark_as_advanced(YAML_INCLUDES YAML_LIBRARIES) diff --git a/glue-codes/openfast-cpp/CMakeLists.txt b/glue-codes/openfast-cpp/CMakeLists.txt index c4dab07b42..762faa6a65 100644 --- a/glue-codes/openfast-cpp/CMakeLists.txt +++ b/glue-codes/openfast-cpp/CMakeLists.txt @@ -25,9 +25,9 @@ find_package(MPI REQUIRED) find_package(LibXml2 REQUIRED) find_package(ZLIB REQUIRED) find_package(HDF5 REQUIRED COMPONENTS C HL) -find_package(YAMLCPP REQUIRED) +find_package(yaml-cpp REQUIRED) -include_directories(${YAML_INCLUDES}) +include_directories(${YAML_CPP_INCLUDE_DIRS}) include_directories(${HDF5_INCLUDES}) include_directories(${HDF5_INCLUDE_DIR}) include_directories(${ZLIB_INCLUDES}) @@ -49,11 +49,10 @@ target_link_libraries(openfastcpplib ${MPI_LIBRARIES} ${CMAKE_DL_LIBS}) -add_executable(openfastcpp - src/FAST_Prog.cpp) - +add_executable(openfastcpp src/FAST_Prog.cpp) target_link_libraries(openfastcpp openfastcpplib openfastlib - ${MPI_LIBRARIES} ${YAML_LIBRARIES} + ${MPI_LIBRARIES} + ${YAML_CPP_LIBRARIES} ${HDF5_C_LIBRARIES} ${HDF5_HL_LIBRARIES} ${ZLIB_LIBRARIES} diff --git a/glue-codes/openfast-cpp/src/FAST_Prog.cpp b/glue-codes/openfast-cpp/src/FAST_Prog.cpp index b7af795629..e6083f7ff3 100644 --- a/glue-codes/openfast-cpp/src/FAST_Prog.cpp +++ b/glue-codes/openfast-cpp/src/FAST_Prog.cpp @@ -3,147 +3,155 @@ #include #include +#include + inline bool checkFileExists(const std::string& name) { - struct stat buffer; - return (stat (name.c_str(), &buffer) == 0); + struct stat buffer; + return (stat (name.c_str(), &buffer) == 0); } void readTurbineData(int iTurb, fast::fastInputs & fi, YAML::Node turbNode) { - - //Read turbine data for a given turbine using the YAML node - fi.globTurbineData[iTurb].TurbID = turbNode["turb_id"].as(); - fi.globTurbineData[iTurb].FASTInputFileName = turbNode["FAST_input_filename"].as() ; - fi.globTurbineData[iTurb].FASTRestartFileName = turbNode["restart_filename"].as() ; - if (turbNode["turbine_base_pos"].IsSequence() ) { - fi.globTurbineData[iTurb].TurbineBasePos = turbNode["turbine_base_pos"].as >() ; - } - if (turbNode["turbine_hub_pos"].IsSequence() ) { - fi.globTurbineData[iTurb].TurbineHubPos = turbNode["turbine_hub_pos"].as >() ; - } - fi.globTurbineData[iTurb].numForcePtsBlade = turbNode["num_force_pts_blade"].as(); - fi.globTurbineData[iTurb].numForcePtsTwr = turbNode["num_force_pts_tower"].as(); - if (turbNode["nacelle_cd"]) {fi.globTurbineData[iTurb].nacelle_cd = turbNode["nacelle_cd"].as();} - if (turbNode["nacelle_area"]) {fi.globTurbineData[iTurb].nacelle_area = turbNode["nacelle_area"].as();} - if (turbNode["air_density"]) {fi.globTurbineData[iTurb].air_density = turbNode["air_density"].as();} + //Read turbine data for a given turbine using the YAML node + fi.globTurbineData[iTurb].TurbID = turbNode["turb_id"].as(); + fi.globTurbineData[iTurb].FASTInputFileName = turbNode["FAST_input_filename"].as(); + fi.globTurbineData[iTurb].FASTRestartFileName = turbNode["restart_filename"].as(); + if (turbNode["turbine_base_pos"].IsSequence() ) { + fi.globTurbineData[iTurb].TurbineBasePos = turbNode["turbine_base_pos"].as >(); + } + if (turbNode["turbine_hub_pos"].IsSequence() ) { + fi.globTurbineData[iTurb].TurbineHubPos = turbNode["turbine_hub_pos"].as >(); + } + fi.globTurbineData[iTurb].numForcePtsBlade = turbNode["num_force_pts_blade"].as(); + fi.globTurbineData[iTurb].numForcePtsTwr = turbNode["num_force_pts_tower"].as(); + if (turbNode["nacelle_cd"]) fi.globTurbineData[iTurb].nacelle_cd = turbNode["nacelle_cd"].as(); + if (turbNode["nacelle_area"]) fi.globTurbineData[iTurb].nacelle_area = turbNode["nacelle_area"].as(); + if (turbNode["air_density"]) fi.globTurbineData[iTurb].air_density = turbNode["air_density"].as(); } void readInputFile(fast::fastInputs & fi, std::string cInterfaceInputFile, double * tEnd) { - fi.comm = MPI_COMM_WORLD; - - // Check if the input file exists and read it - if ( checkFileExists(cInterfaceInputFile) ) { - - YAML::Node cDriverInp = YAML::LoadFile(cInterfaceInputFile); - - fi.nTurbinesGlob = cDriverInp["nTurbinesGlob"].as(); - - if (fi.nTurbinesGlob > 0) { - - if(cDriverInp["dryRun"]) { - fi.dryRun = cDriverInp["dryRun"].as(); - } - - 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(); - fi.nEveryCheckPoint = cDriverInp["nEveryCheckPoint"].as(); - fi.dtFAST = cDriverInp["dtFAST"].as(); - fi.tMax = cDriverInp["tMax"].as(); // tMax is the total duration to which you want to run FAST. This should be the same or greater than the max time given in the FAST fst file. Choose this carefully as FAST writes the output file only at this point if you choose the binary file output. - - if(cDriverInp["superController"]) { - fi.scStatus = cDriverInp["superController"].as(); - fi.scLibFile = cDriverInp["scLibFile"].as(); - fi.numScInputs = cDriverInp["numScInputs"].as(); - fi.numScOutputs = cDriverInp["numScOutputs"].as(); - } - - fi.globTurbineData.resize(fi.nTurbinesGlob); - for (int iTurb=0; iTurb < fi.nTurbinesGlob; iTurb++) { - if (cDriverInp["Turbine" + std::to_string(iTurb)]) { - readTurbineData(iTurb, fi, cDriverInp["Turbine" + std::to_string(iTurb)] ); - } else { - throw std::runtime_error("Node for Turbine" + std::to_string(iTurb) + " not present in input file or I cannot read it"); - } - } - + fi.comm = MPI_COMM_WORLD; + + // Check if the input file exists and read it + if ( checkFileExists(cInterfaceInputFile) ) { + + YAML::Node cDriverInp = YAML::LoadFile(cInterfaceInputFile); + + fi.nTurbinesGlob = cDriverInp["nTurbinesGlob"].as(); + + if (fi.nTurbinesGlob > 0) { + + if(cDriverInp["dryRun"]) { + fi.dryRun = cDriverInp["dryRun"].as(); + } + + 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(); + fi.nEveryCheckPoint = cDriverInp["nEveryCheckPoint"].as(); + fi.dtFAST = cDriverInp["dtFAST"].as(); + fi.tMax = cDriverInp["tMax"].as(); // tMax is the total duration to which you want to run FAST. This should be the same or greater than the max time given in the FAST fst file. Choose this carefully as FAST writes the output file only at this point if you choose the binary file output. + + if(cDriverInp["superController"]) { + fi.scStatus = cDriverInp["superController"].as(); + fi.scLibFile = cDriverInp["scLibFile"].as(); + fi.numScInputs = cDriverInp["numScInputs"].as(); + fi.numScOutputs = cDriverInp["numScOutputs"].as(); + } + + fi.globTurbineData.resize(fi.nTurbinesGlob); + for (int iTurb=0; iTurb < fi.nTurbinesGlob; iTurb++) { + if (cDriverInp["Turbine" + std::to_string(iTurb)]) { + readTurbineData(iTurb, fi, cDriverInp["Turbine" + std::to_string(iTurb)] ); + } else { + throw std::runtime_error("Node for Turbine" + std::to_string(iTurb) + " not present in input file or I cannot read it"); + } + } + + } else { + throw std::runtime_error("Number of turbines <= 0 "); + } + } else { - throw std::runtime_error("Number of turbines <= 0 "); + throw std::runtime_error("Input file " + cInterfaceInputFile + " does not exist or I cannot access it"); } - - } else { - throw std::runtime_error("Input file " + cInterfaceInputFile + " does not exist or I cannot access it"); - } - } -int main() { - int iErr; - int nProcs; - int rank; - std::vector torque (3, 0.0); - std::vector thrust (3, 0.0); - - iErr = MPI_Init(NULL, NULL); - iErr = MPI_Comm_size( MPI_COMM_WORLD, &nProcs); - iErr = MPI_Comm_rank( MPI_COMM_WORLD, &rank); - - double tEnd ; // This doesn't belong in the FAST - C++ interface - int ntEnd ; // This doesn't belong in the FAST - C++ interface - - std::string cDriverInputFile="cDriver.i"; - fast::OpenFAST FAST; - fast::fastInputs fi ; - try { - readInputFile(fi, cDriverInputFile, &tEnd); - } - catch( const std::runtime_error & ex) { - std::cerr << ex.what() << std::endl ; - std::cerr << "Program quitting now" << std::endl ; - return 1; - } - ntEnd = tEnd/fi.dtFAST; //Calculate the last time step - - FAST.setInputs(fi); - FAST.allocateTurbinesToProcsSimple(); - // Or allocate turbines to procs by calling "setTurbineProcNo(iTurbGlob, procId)" for turbine. - - FAST.init(); - if (FAST.isTimeZero()) { - FAST.solution0(); - } - - if( !FAST.isDryRun() ) { +int main(int argc, char** argv) { + if (argc != 2) { + std::cerr << "Incorrect syntax. Try: openfastcpp inputfile.yaml" << std::endl ; + return 1; + } + + int iErr; + int nProcs; + int rank; + std::vector torque (3, 0.0); + std::vector thrust (3, 0.0); + + iErr = MPI_Init(NULL, NULL); + iErr = MPI_Comm_size( MPI_COMM_WORLD, &nProcs); + iErr = MPI_Comm_rank( MPI_COMM_WORLD, &rank); + + double tEnd ; // This doesn't belong in the FAST - C++ interface + int ntEnd ; // This doesn't belong in the FAST - C++ interface + + std::string cDriverInputFile=argv[1]; + fast::OpenFAST FAST; + fast::fastInputs fi ; + try { + readInputFile(fi, cDriverInputFile, &tEnd); + } catch( const std::runtime_error & ex) { + std::cerr << ex.what() << std::endl ; + std::cerr << "Program quitting now" << std::endl ; + return 1; + } + + // Calculate the last time step + ntEnd = tEnd/fi.dtFAST; + + FAST.setInputs(fi); + FAST.allocateTurbinesToProcsSimple(); + // Or allocate turbines to procs by calling "setTurbineProcNo(iTurbGlob, procId)" for turbine. + + FAST.init(); + + if (FAST.isTimeZero()) FAST.solution0(); + + if ( FAST.isDryRun() ) { + FAST.end() ; + MPI_Finalize() ; + return 0; + } + for (int nt = FAST.get_ntStart(); nt < ntEnd; nt++) { - FAST.step(); - if (FAST.isDebug()) { - FAST.computeTorqueThrust(0,torque,thrust); - std::cout.precision(16); - std::cout << "Torque = " << torque[0] << " " << torque[1] << " " << torque[2] << std::endl ; - std::cout << "Thrust = " << thrust[0] << " " << thrust[1] << " " << thrust[2] << std::endl ; - } + FAST.step(); + if (FAST.isDebug()) { + FAST.computeTorqueThrust(0,torque,thrust); + std::cout.precision(16); + std::cout << "Torque = " << torque[0] << " " << torque[1] << " " << torque[2] << std::endl ; + std::cout << "Thrust = " << thrust[0] << " " << thrust[1] << " " << thrust[2] << std::endl ; + } } - } - FAST.end() ; - MPI_Finalize() ; + FAST.end() ; + MPI_Finalize() ; - return 0; - -} + return 0; +} diff --git a/glue-codes/openfast-cpp/src/OpenFAST.cpp b/glue-codes/openfast-cpp/src/OpenFAST.cpp index c2901fec4f..b4c220d7d8 100644 --- a/glue-codes/openfast-cpp/src/OpenFAST.cpp +++ b/glue-codes/openfast-cpp/src/OpenFAST.cpp @@ -24,7 +24,6 @@ numScOutputs(0) //Nothing to do here } - //Constructor fast::OpenFAST::OpenFAST(): nTurbinesGlob(0), @@ -35,638 +34,622 @@ timeZero(false) { } +fast::OpenFAST::~OpenFAST(){ } + inline bool fast::OpenFAST::checkFileExists(const std::string& name) { - struct stat buffer; - return (stat (name.c_str(), &buffer) == 0); + struct stat buffer; + return (stat (name.c_str(), &buffer) == 0); } void fast::OpenFAST::init() { - // Temporary buffer to pass filenames to OpenFAST fortran subroutines - char currentFileName[INTERFACE_STRING_LENGTH]; - - allocateMemory(); - - if (!dryRun) { - switch (simStart) { - - case fast::trueRestart: + // Temporary buffer to pass filenames to OpenFAST fortran subroutines + char currentFileName[INTERFACE_STRING_LENGTH]; + + allocateMemory(); + + if (!dryRun) { + switch (simStart) { + + case fast::trueRestart: + + for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { + /* note that this will set nt_global inside the FAST library */ + std::copy( + CheckpointFileRoot[iTurb].data(), + CheckpointFileRoot[iTurb].data() + (CheckpointFileRoot[iTurb].size() + 1), + currentFileName + ); + FAST_OpFM_Restart( + &iTurb, + currentFileName, + &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 (nTurbinesProc > 0) velNodeDataFile = openVelocityDataFile(false); + + if(scStatus) { + sc->readRestartFile(nt_global); + } + + break ; + + case fast::init: + + // this calls the Init() routines of each module + + for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { + std::copy( + FASTInputFileName[iTurb].data(), + FASTInputFileName[iTurb].data() + (FASTInputFileName[iTurb].size() + 1), + currentFileName + ); + FAST_OpFM_Init( + &iTurb, + &tMax, + currentFileName, + &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; + if(numVelPtsTwr[iTurb] == 0) { + numForcePtsTwr[iTurb] = 0; + std::cout << "Aerodyn doesn't want to calculate forces on the tower. All actuator points on the tower are turned off for turbine " << turbineMapProcToGlob[iTurb] << "." << std::endl ; + } + + 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++) { + std::copy( + FASTInputFileName[iTurb].data(), + FASTInputFileName[iTurb].data() + (FASTInputFileName[iTurb].size() + 1), + currentFileName + ); + FAST_OpFM_Init( + &iTurb, + &tMax, + currentFileName, + &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 ; + } + } + } + + int nTimesteps; + + if (nTurbinesProc > 0) { + readVelocityData(ntStart); + } + for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { + applyVelocityData(0, iTurb, cDriver_Output_to_FAST[iTurb], velNodeData[iTurb]); + } + solution0() ; - for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { - /* note that this will set nt_global inside the FAST library */ - std::copy(CheckpointFileRoot[iTurb].data(), - CheckpointFileRoot[iTurb].data() + (CheckpointFileRoot[iTurb].size() + 1), - currentFileName); - FAST_OpFM_Restart( - &iTurb, currentFileName, &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; + 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(); + } - 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); - } + break; - if (nTurbinesProc > 0) velNodeDataFile = openVelocityDataFile(false); + case fast::simStartType_END: - if(scStatus) { - sc->readRestartFile(nt_global); - } - - break ; - - case fast::init: - - // this calls the Init() routines of each module - - for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { - std::copy(FASTInputFileName[iTurb].data(), - FASTInputFileName[iTurb].data() + (FASTInputFileName[iTurb].size() + 1), - currentFileName); - FAST_OpFM_Init(&iTurb, &tMax, currentFileName, &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; - if(numVelPtsTwr[iTurb] == 0) { - numForcePtsTwr[iTurb] = 0; - std::cout << "Aerodyn doesn't want to calculate forces on the tower. All actuator points on the tower are turned off for turbine " << turbineMapProcToGlob[iTurb] << "." << std::endl ; - } - - - 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++) { - std::copy(FASTInputFileName[iTurb].data(), - FASTInputFileName[iTurb].data() + (FASTInputFileName[iTurb].size() + 1), - currentFileName); - FAST_OpFM_Init(&iTurb, &tMax, currentFileName, &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 ; - } - } - } - - int nTimesteps; - - if (nTurbinesProc > 0) { - readVelocityData(ntStart); - } - 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; + break; + } } - - } } 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]); - // } - - if(scStatus) { - - sc->init(nTurbinesGlob, numScInputs, numScOutputs); - - sc->calcOutputs(scOutputsGlob); - fillScOutputsLoc(); - } + 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++) { + if(scStatus) { - FAST_OpFM_Solution0(&iTurb, &ErrStat, ErrMsg); - checkError(ErrStat, ErrMsg); + sc->init(nTurbinesGlob, numScInputs, numScOutputs); - } + sc->calcOutputs(scOutputsGlob); + fillScOutputsLoc(); + } - timeZero = false; + for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { + FAST_OpFM_Solution0(&iTurb, &ErrStat, ErrMsg); + checkError(ErrStat, ErrMsg); + } - if (scStatus) { - fillScInputsGlob(); // Update inputs to super controller - } - } + timeZero = false; + if (scStatus) { + fillScInputsGlob(); // Update inputs to super controller + } + } } void fast::OpenFAST::step() { - /* ****************************** + /* ****************************** 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) - - writeVelocityData(velNodeDataFile, iTurb, nt_global, cDriver_Input_from_FAST[iTurb], cDriver_Output_to_FAST[iTurb]); - - if ( isDebug() ) { - - std::ofstream fastcpp_velocity_file; - fastcpp_velocity_file.open("fastcpp_velocity.csv") ; - 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].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); - - // Compute the force from the nacelle only if the drag coefficient is - // greater than zero - if (nacelle_cd[iTurb]>0.) { - - calc_nacelle_force ( - - cDriver_Output_to_FAST[iTurb].u[0], - cDriver_Output_to_FAST[iTurb].v[0], - cDriver_Output_to_FAST[iTurb].w[0], - nacelle_cd[iTurb], - nacelle_area[iTurb], - air_density[iTurb], - cDriver_Input_from_FAST[iTurb].fx[0], - cDriver_Input_from_FAST[iTurb].fy[0], - cDriver_Input_from_FAST[iTurb].fz[0] - - ); - - } - - if ( isDebug() ) { - std::ofstream actuatorForcesFile; - 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 ; - } - actuatorForcesFile.close() ; - } - - } - - 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; - - if ( (((nt_global - ntStart) % nEveryCheckPoint) == 0 ) && (nt_global != ntStart) ) { - // Use default FAST naming convention for checkpoint file - // . - char dummyCheckPointRoot[INTERFACE_STRING_LENGTH] = " "; - // Ensure that we have a null character - dummyCheckPointRoot[1] = 0; + ********************************* */ - if (nTurbinesProc > 0) backupVelocityDataFile(nt_global, velNodeDataFile); + if(scStatus) { + sc->calcOutputs(scOutputsGlob); + fillScOutputsLoc(); + } for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { - FAST_CreateCheckpoint(&iTurb, dummyCheckPointRoot, &ErrStat, ErrMsg); - checkError(ErrStat, ErrMsg); + + // 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; + fastcpp_velocity_file.open("fastcpp_velocity.csv") ; + 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].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); + + // Compute the force from the nacelle only if the drag coefficient is + // greater than zero + if (nacelle_cd[iTurb]>0.) { + calc_nacelle_force ( + cDriver_Output_to_FAST[iTurb].u[0], + cDriver_Output_to_FAST[iTurb].v[0], + cDriver_Output_to_FAST[iTurb].w[0], + nacelle_cd[iTurb], + nacelle_area[iTurb], + air_density[iTurb], + cDriver_Input_from_FAST[iTurb].fx[0], + cDriver_Input_from_FAST[iTurb].fy[0], + cDriver_Input_from_FAST[iTurb].fz[0] + ); + } + + if ( isDebug() ) { + std::ofstream actuatorForcesFile; + 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 ; + } + actuatorForcesFile.close() ; + } } + if(scStatus) { - if (fastMPIRank == 0) { - sc->writeRestartFile(nt_global); - } + 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; + + if ( (((nt_global - ntStart) % nEveryCheckPoint) == 0 ) && (nt_global != ntStart) ) { + // Use default FAST naming convention for checkpoint file + // . + char dummyCheckPointRoot[INTERFACE_STRING_LENGTH] = " "; + // Ensure that we have a null character + dummyCheckPointRoot[1] = 0; + + if (nTurbinesProc > 0) backupVelocityDataFile(nt_global, velNodeDataFile); + + for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { + FAST_CreateCheckpoint(&iTurb, dummyCheckPointRoot, &ErrStat, ErrMsg); + checkError(ErrStat, ErrMsg); + } + if(scStatus) { + if (fastMPIRank == 0) { + sc->writeRestartFile(nt_global); + } + } + } } void fast::OpenFAST::stepNoWrite() { - /* ****************************** - set inputs from this code and call FAST: - ********************************* */ + /* ****************************** + set inputs from this code and call FAST: + ********************************* */ - if(scStatus) { - sc->calcOutputs(scOutputsGlob); - fillScOutputsLoc(); - } + if(scStatus) { + sc->calcOutputs(scOutputsGlob); + fillScOutputsLoc(); + } - for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { + for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { - // set wind speeds at original locations - // setOutputsToFAST(cDriver_Input_from_FAST[iTurb], cDriver_Output_to_FAST[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); + // 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' - } + 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; - + nt_global = nt_global + 1; } -fast::OpenFAST::~OpenFAST(){ -} +void fast::OpenFAST::calc_nacelle_force(const float & u, const float & v, const float & w, const float & cd, const float & area, const float & rho, float & fx, float & fy, float & fz) { + // Calculate the force on the nacelle (fx,fy,fz) given the + // velocity sampled at the nacelle point (u,v,w), + // drag coefficient 'cd' and nacelle area 'area' -void fast::OpenFAST::calc_nacelle_force( - const float & u, - const float & v, - const float & w, - const float & cd, - const float & area, - const float & rho, - float & fx, - float & fy, - float & fz) { - // Calculate the force on the nacelle (fx,fy,fz) given the - // velocity sampled at the nacelle point (u,v,w), - // drag coefficient 'cd' and nacelle area 'area' - - // The velocity magnitude - float Vmag = std::sqrt(u * u + v * v + w * w); - - // Velocity correction based on Martinez-Tossas PhD Thesis 2017 - // The correction samples the velocity at the center of the - // Gaussian kernel and scales it to obtain the inflow velocity - float epsilon_d = std::sqrt(2.0 / M_PI * cd * area); - float correction = 1. / (1.0 - cd * area / - (4.0 * M_PI * epsilon_d * epsilon_d)); - - // Compute the force for each velocity component - fx = rho * 1./2. * cd * area * Vmag * u * correction * correction; - fy = rho * 1./2. * cd * area * Vmag * v * correction * correction; - fz = rho * 1./2. * cd * area * Vmag * w * correction * correction; - } + // The velocity magnitude + float Vmag = std::sqrt(u * u + v * v + w * w); -void fast::OpenFAST::setInputs(const fast::fastInputs & fi ) { + // Velocity correction based on Martinez-Tossas PhD Thesis 2017 + // The correction samples the velocity at the center of the + // Gaussian kernel and scales it to obtain the inflow velocity + float epsilon_d = std::sqrt(2.0 / M_PI * cd * area); + float correction = 1. / (1.0 - cd * area / (4.0 * M_PI * epsilon_d * epsilon_d)); + + // Compute the force for each velocity component + fx = rho * 1./2. * cd * area * Vmag * u * correction * correction; + fy = rho * 1./2. * cd * area * Vmag * v * correction * correction; + fz = rho * 1./2. * cd * area * Vmag * w * correction * correction; +} +void fast::OpenFAST::setInputs(const fast::fastInputs & fi ) { - mpiComm = fi.comm; + mpiComm = fi.comm; - MPI_Comm_rank(mpiComm, &worldMPIRank); - MPI_Comm_group(mpiComm, &worldMPIGroup); + MPI_Comm_rank(mpiComm, &worldMPIRank); + MPI_Comm_group(mpiComm, &worldMPIGroup); - nTurbinesGlob = fi.nTurbinesGlob; + nTurbinesGlob = fi.nTurbinesGlob; if (nTurbinesGlob > 0) { - - dryRun = fi.dryRun; - - debug = fi.debug; - - tStart = fi.tStart; - simStart = fi.simStart; - nEveryCheckPoint = fi.nEveryCheckPoint; - tMax = fi.tMax; - loadSuperController(fi); - dtFAST = fi.dtFAST; - - ntStart = int(tStart/dtFAST); - - if (simStart == fast::restartDriverInitFAST) { - nt_global = 0; - } else { - nt_global = ntStart; - } - - globTurbineData.resize(nTurbinesGlob); - globTurbineData = fi.globTurbineData; - } else { - throw std::runtime_error("Number of turbines < 0 "); - } - -} + dryRun = fi.dryRun; + debug = fi.debug; -void fast::OpenFAST::checkError(const int ErrStat, const char * ErrMsg){ + tStart = fi.tStart; + simStart = fi.simStart; + nEveryCheckPoint = fi.nEveryCheckPoint; + tMax = fi.tMax; + loadSuperController(fi); + dtFAST = fi.dtFAST; - if (ErrStat != ErrID_None){ + ntStart = int(tStart/dtFAST); - if (ErrStat >= AbortErrLev){ - throw std::runtime_error(ErrMsg); - } + if (simStart == fast::restartDriverInitFAST) { + nt_global = 0; + } else { + nt_global = ntStart; + } - } + globTurbineData.resize(nTurbinesGlob); + globTurbineData = fi.globTurbineData; + } else { + throw std::runtime_error("Number of turbines < 0 "); + } } -void fast::OpenFAST::setOutputsToFAST(OpFM_InputType_t cDriver_Input_from_FAST, OpFM_OutputType_t cDriver_Output_to_FAST){ +void fast::OpenFAST::checkError(const int ErrStat, const char * ErrMsg){ + if (ErrStat != ErrID_None){ + if (ErrStat >= AbortErrLev){ + throw std::runtime_error(ErrMsg); + } + } +} - // routine sets the u-v-w wind speeds used in FAST and the SuperController inputs +void fast::OpenFAST::setOutputsToFAST(OpFM_InputType_t cDriver_Input_from_FAST, OpFM_OutputType_t cDriver_Output_to_FAST){ - for (int j = 0; j < cDriver_Output_to_FAST.u_Len; j++){ - cDriver_Output_to_FAST.u[j] = (float) 10.0*pow((cDriver_Input_from_FAST.pzVel[j] / 90.0), 0.2); // 0.2 power law wind profile using reference 10 m/s at 90 meters - cDriver_Output_to_FAST.v[j] = 0.0; - cDriver_Output_to_FAST.w[j] = 0.0; - } + // routine sets the u-v-w wind speeds used in FAST and the SuperController inputs - // // 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) - // } + for (int j = 0; j < cDriver_Output_to_FAST.u_Len; j++){ + cDriver_Output_to_FAST.u[j] = (float) 10.0*pow((cDriver_Input_from_FAST.pzVel[j] / 90.0), 0.2); // 0.2 power law wind profile using reference 10 m/s at 90 meters + cDriver_Output_to_FAST.v[j] = 0.0; + 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) + // } } void fast::OpenFAST::getApproxHubPos(double* currentCoords, int iTurbGlob, int nSize) { - - assert(nSize==3); - // Get hub position of Turbine 'iTurbGlob' - for(int i =0; i rDistForce(nForcePtsBlade) ; - for(int j=0; j < nForcePtsBlade; j++) { - int iNodeForce = 1 + iBlade * nForcePtsBlade + j ; //The number of actuator force points is always the same for all blades - rDistForce[j] = std::sqrt( - (cDriver_Input_from_FAST[iTurb].pxForce[iNodeForce] - cDriver_Input_from_FAST[iTurb].pxForce[0])*(cDriver_Input_from_FAST[iTurb].pxForce[iNodeForce] - cDriver_Input_from_FAST[iTurb].pxForce[0]) - + (cDriver_Input_from_FAST[iTurb].pyForce[iNodeForce] - cDriver_Input_from_FAST[iTurb].pyForce[0])*(cDriver_Input_from_FAST[iTurb].pyForce[iNodeForce] - cDriver_Input_from_FAST[iTurb].pyForce[0]) - + (cDriver_Input_from_FAST[iTurb].pzForce[iNodeForce] - cDriver_Input_from_FAST[iTurb].pzForce[0])*(cDriver_Input_from_FAST[iTurb].pzForce[iNodeForce] - cDriver_Input_from_FAST[iTurb].pzForce[0]) - ); - } - - // Interpolate to the velocity nodes - int nVelPtsBlade = get_numVelPtsBladeLoc(iTurb); - for(int j=0; j < nVelPtsBlade; j++) { - int iNodeVel = 1 + iBlade * nVelPtsBlade + j ; //Assumes the same number of velocity (Aerodyn) nodes for all blades - double rDistVel = std::sqrt( - (cDriver_Input_from_FAST[iTurb].pxVel[iNodeVel] - cDriver_Input_from_FAST[iTurb].pxVel[0])*(cDriver_Input_from_FAST[iTurb].pxVel[iNodeVel] - cDriver_Input_from_FAST[iTurb].pxVel[0]) - + (cDriver_Input_from_FAST[iTurb].pyVel[iNodeVel] - cDriver_Input_from_FAST[iTurb].pyVel[0])*(cDriver_Input_from_FAST[iTurb].pyVel[iNodeVel] - cDriver_Input_from_FAST[iTurb].pyVel[0]) - + (cDriver_Input_from_FAST[iTurb].pzVel[iNodeVel] - cDriver_Input_from_FAST[iTurb].pzVel[0])*(cDriver_Input_from_FAST[iTurb].pzVel[iNodeVel] - cDriver_Input_from_FAST[iTurb].pzVel[0]) - ); - //Find nearest two force nodes - int jForceLower = 0; - while ( (rDistForce[jForceLower+1] < rDistVel) && ( jForceLower < (nForcePtsBlade-2)) ) { - jForceLower = jForceLower + 1; - } - int iNodeForceLower = 1 + iBlade * nForcePtsBlade + jForceLower ; - double rInterp = (rDistVel - rDistForce[jForceLower])/(rDistForce[jForceLower+1]-rDistForce[jForceLower]); - cDriver_Output_to_FAST[iTurb].u[iNodeVel] = forceNodeVel[iTurb][iNodeForceLower][0] + rInterp * (forceNodeVel[iTurb][iNodeForceLower+1][0] - forceNodeVel[iTurb][iNodeForceLower][0] ); - cDriver_Output_to_FAST[iTurb].v[iNodeVel] = forceNodeVel[iTurb][iNodeForceLower][1] + rInterp * (forceNodeVel[iTurb][iNodeForceLower+1][1] - forceNodeVel[iTurb][iNodeForceLower][1] ); - cDriver_Output_to_FAST[iTurb].w[iNodeVel] = forceNodeVel[iTurb][iNodeForceLower][2] + rInterp * (forceNodeVel[iTurb][iNodeForceLower+1][2] - forceNodeVel[iTurb][iNodeForceLower][2] ); - } - } + // Do the blades first + int nBlades = get_numBladesLoc(iTurb); + for(int iBlade=0; iBlade < nBlades; iBlade++) { + // Create interpolating parameter - Distance from hub + int nForcePtsBlade = get_numForcePtsBladeLoc(iTurb); + std::vector rDistForce(nForcePtsBlade) ; + for(int j=0; j < nForcePtsBlade; j++) { + int iNodeForce = 1 + iBlade * nForcePtsBlade + j ; //The number of actuator force points is always the same for all blades + rDistForce[j] = std::sqrt( + (cDriver_Input_from_FAST[iTurb].pxForce[iNodeForce] - cDriver_Input_from_FAST[iTurb].pxForce[0])*(cDriver_Input_from_FAST[iTurb].pxForce[iNodeForce] - cDriver_Input_from_FAST[iTurb].pxForce[0]) + + (cDriver_Input_from_FAST[iTurb].pyForce[iNodeForce] - cDriver_Input_from_FAST[iTurb].pyForce[0])*(cDriver_Input_from_FAST[iTurb].pyForce[iNodeForce] - cDriver_Input_from_FAST[iTurb].pyForce[0]) + + (cDriver_Input_from_FAST[iTurb].pzForce[iNodeForce] - cDriver_Input_from_FAST[iTurb].pzForce[0])*(cDriver_Input_from_FAST[iTurb].pzForce[iNodeForce] - cDriver_Input_from_FAST[iTurb].pzForce[0]) + ); + } + + // Interpolate to the velocity nodes + int nVelPtsBlade = get_numVelPtsBladeLoc(iTurb); + for(int j=0; j < nVelPtsBlade; j++) { + int iNodeVel = 1 + iBlade * nVelPtsBlade + j ; //Assumes the same number of velocity (Aerodyn) nodes for all blades + double rDistVel = std::sqrt( + (cDriver_Input_from_FAST[iTurb].pxVel[iNodeVel] - cDriver_Input_from_FAST[iTurb].pxVel[0])*(cDriver_Input_from_FAST[iTurb].pxVel[iNodeVel] - cDriver_Input_from_FAST[iTurb].pxVel[0]) + + (cDriver_Input_from_FAST[iTurb].pyVel[iNodeVel] - cDriver_Input_from_FAST[iTurb].pyVel[0])*(cDriver_Input_from_FAST[iTurb].pyVel[iNodeVel] - cDriver_Input_from_FAST[iTurb].pyVel[0]) + + (cDriver_Input_from_FAST[iTurb].pzVel[iNodeVel] - cDriver_Input_from_FAST[iTurb].pzVel[0])*(cDriver_Input_from_FAST[iTurb].pzVel[iNodeVel] - cDriver_Input_from_FAST[iTurb].pzVel[0]) + ); + //Find nearest two force nodes + int jForceLower = 0; + while ( (rDistForce[jForceLower+1] < rDistVel) && ( jForceLower < (nForcePtsBlade-2)) ) { + jForceLower = jForceLower + 1; + } + int iNodeForceLower = 1 + iBlade * nForcePtsBlade + jForceLower ; + double rInterp = (rDistVel - rDistForce[jForceLower])/(rDistForce[jForceLower+1]-rDistForce[jForceLower]); + cDriver_Output_to_FAST[iTurb].u[iNodeVel] = forceNodeVel[iTurb][iNodeForceLower][0] + rInterp * (forceNodeVel[iTurb][iNodeForceLower+1][0] - forceNodeVel[iTurb][iNodeForceLower][0] ); + cDriver_Output_to_FAST[iTurb].v[iNodeVel] = forceNodeVel[iTurb][iNodeForceLower][1] + rInterp * (forceNodeVel[iTurb][iNodeForceLower+1][1] - forceNodeVel[iTurb][iNodeForceLower][1] ); + cDriver_Output_to_FAST[iTurb].w[iNodeVel] = forceNodeVel[iTurb][iNodeForceLower][2] + rInterp * (forceNodeVel[iTurb][iNodeForceLower+1][2] - forceNodeVel[iTurb][iNodeForceLower][2] ); + } + } - // Now the tower if present and used - int nVelPtsTower = get_numVelPtsTwrLoc(iTurb); - if ( nVelPtsTower > 0 ) { - - // Create interpolating parameter - Distance from first node from ground - int nForcePtsTower = get_numForcePtsTwrLoc(iTurb); - std::vector hDistForce(nForcePtsTower) ; - int iNodeBotTowerForce = 1 + nBlades * get_numForcePtsBladeLoc(iTurb); // The number of actuator force points is always the same for all blades - for(int j=0; j < nForcePtsTower; j++) { - int iNodeForce = iNodeBotTowerForce + j ; - hDistForce[j] = std::sqrt( - (cDriver_Input_from_FAST[iTurb].pxForce[iNodeForce] - cDriver_Input_from_FAST[iTurb].pxForce[iNodeBotTowerForce])*(cDriver_Input_from_FAST[iTurb].pxForce[iNodeForce] - cDriver_Input_from_FAST[iTurb].pxForce[iNodeBotTowerForce]) - + (cDriver_Input_from_FAST[iTurb].pyForce[iNodeForce] - cDriver_Input_from_FAST[iTurb].pyForce[iNodeBotTowerForce])*(cDriver_Input_from_FAST[iTurb].pyForce[iNodeForce] - cDriver_Input_from_FAST[iTurb].pyForce[iNodeBotTowerForce]) - + (cDriver_Input_from_FAST[iTurb].pzForce[iNodeForce] - cDriver_Input_from_FAST[iTurb].pzForce[iNodeBotTowerForce])*(cDriver_Input_from_FAST[iTurb].pzForce[iNodeForce] - cDriver_Input_from_FAST[iTurb].pzForce[iNodeBotTowerForce]) - ); - } - - - int iNodeBotTowerVel = 1 + nBlades * get_numVelPtsBladeLoc(iTurb); // Assumes the same number of velocity (Aerodyn) nodes for all blades - for(int j=0; j < nVelPtsTower; j++) { - int iNodeVel = iNodeBotTowerVel + j ; - double hDistVel = std::sqrt( - (cDriver_Input_from_FAST[iTurb].pxVel[iNodeVel] - cDriver_Input_from_FAST[iTurb].pxVel[iNodeBotTowerVel])*(cDriver_Input_from_FAST[iTurb].pxVel[iNodeVel] - cDriver_Input_from_FAST[iTurb].pxVel[iNodeBotTowerVel]) - + (cDriver_Input_from_FAST[iTurb].pyVel[iNodeVel] - cDriver_Input_from_FAST[iTurb].pyVel[iNodeBotTowerVel])*(cDriver_Input_from_FAST[iTurb].pyVel[iNodeVel] - cDriver_Input_from_FAST[iTurb].pyVel[iNodeBotTowerVel]) - + (cDriver_Input_from_FAST[iTurb].pzVel[iNodeVel] - cDriver_Input_from_FAST[iTurb].pzVel[iNodeBotTowerVel])*(cDriver_Input_from_FAST[iTurb].pzVel[iNodeVel] - cDriver_Input_from_FAST[iTurb].pzVel[iNodeBotTowerVel]) - ); - //Find nearest two force nodes - int jForceLower = 0; - while ( (hDistForce[jForceLower+1] < hDistVel) && ( jForceLower < (nForcePtsTower-2)) ) { - jForceLower = jForceLower + 1; - } - int iNodeForceLower = iNodeBotTowerForce + jForceLower ; - double rInterp = (hDistVel - hDistForce[jForceLower])/(hDistForce[jForceLower+1]-hDistForce[jForceLower]); - cDriver_Output_to_FAST[iTurb].u[iNodeVel] = forceNodeVel[iTurb][iNodeForceLower][0] + rInterp * (forceNodeVel[iTurb][iNodeForceLower+1][0] - forceNodeVel[iTurb][iNodeForceLower][0] ); - cDriver_Output_to_FAST[iTurb].v[iNodeVel] = forceNodeVel[iTurb][iNodeForceLower][1] + rInterp * (forceNodeVel[iTurb][iNodeForceLower+1][1] - forceNodeVel[iTurb][iNodeForceLower][1] ); - cDriver_Output_to_FAST[iTurb].w[iNodeVel] = forceNodeVel[iTurb][iNodeForceLower][2] + rInterp * (forceNodeVel[iTurb][iNodeForceLower+1][2] - forceNodeVel[iTurb][iNodeForceLower][2] ); - } - } - - } - + // Now the tower if present and used + int nVelPtsTower = get_numVelPtsTwrLoc(iTurb); + if ( nVelPtsTower > 0 ) { + + // Create interpolating parameter - Distance from first node from ground + int nForcePtsTower = get_numForcePtsTwrLoc(iTurb); + std::vector hDistForce(nForcePtsTower) ; + int iNodeBotTowerForce = 1 + nBlades * get_numForcePtsBladeLoc(iTurb); // The number of actuator force points is always the same for all blades + for(int j=0; j < nForcePtsTower; j++) { + int iNodeForce = iNodeBotTowerForce + j ; + hDistForce[j] = std::sqrt( + (cDriver_Input_from_FAST[iTurb].pxForce[iNodeForce] - cDriver_Input_from_FAST[iTurb].pxForce[iNodeBotTowerForce])*(cDriver_Input_from_FAST[iTurb].pxForce[iNodeForce] - cDriver_Input_from_FAST[iTurb].pxForce[iNodeBotTowerForce]) + + (cDriver_Input_from_FAST[iTurb].pyForce[iNodeForce] - cDriver_Input_from_FAST[iTurb].pyForce[iNodeBotTowerForce])*(cDriver_Input_from_FAST[iTurb].pyForce[iNodeForce] - cDriver_Input_from_FAST[iTurb].pyForce[iNodeBotTowerForce]) + + (cDriver_Input_from_FAST[iTurb].pzForce[iNodeForce] - cDriver_Input_from_FAST[iTurb].pzForce[iNodeBotTowerForce])*(cDriver_Input_from_FAST[iTurb].pzForce[iNodeForce] - cDriver_Input_from_FAST[iTurb].pzForce[iNodeBotTowerForce]) + ); + } + + int iNodeBotTowerVel = 1 + nBlades * get_numVelPtsBladeLoc(iTurb); // Assumes the same number of velocity (Aerodyn) nodes for all blades + for(int j=0; j < nVelPtsTower; j++) { + int iNodeVel = iNodeBotTowerVel + j ; + double hDistVel = std::sqrt( + (cDriver_Input_from_FAST[iTurb].pxVel[iNodeVel] - cDriver_Input_from_FAST[iTurb].pxVel[iNodeBotTowerVel])*(cDriver_Input_from_FAST[iTurb].pxVel[iNodeVel] - cDriver_Input_from_FAST[iTurb].pxVel[iNodeBotTowerVel]) + + (cDriver_Input_from_FAST[iTurb].pyVel[iNodeVel] - cDriver_Input_from_FAST[iTurb].pyVel[iNodeBotTowerVel])*(cDriver_Input_from_FAST[iTurb].pyVel[iNodeVel] - cDriver_Input_from_FAST[iTurb].pyVel[iNodeBotTowerVel]) + + (cDriver_Input_from_FAST[iTurb].pzVel[iNodeVel] - cDriver_Input_from_FAST[iTurb].pzVel[iNodeBotTowerVel])*(cDriver_Input_from_FAST[iTurb].pzVel[iNodeVel] - cDriver_Input_from_FAST[iTurb].pzVel[iNodeBotTowerVel]) + ); + //Find nearest two force nodes + int jForceLower = 0; + while ( (hDistForce[jForceLower+1] < hDistVel) && ( jForceLower < (nForcePtsTower-2)) ) { + jForceLower = jForceLower + 1; + } + int iNodeForceLower = iNodeBotTowerForce + jForceLower ; + double rInterp = (hDistVel - hDistForce[jForceLower])/(hDistForce[jForceLower+1]-hDistForce[jForceLower]); + cDriver_Output_to_FAST[iTurb].u[iNodeVel] = forceNodeVel[iTurb][iNodeForceLower][0] + rInterp * (forceNodeVel[iTurb][iNodeForceLower+1][0] - forceNodeVel[iTurb][iNodeForceLower][0] ); + cDriver_Output_to_FAST[iTurb].v[iNodeVel] = forceNodeVel[iTurb][iNodeForceLower][1] + rInterp * (forceNodeVel[iTurb][iNodeForceLower+1][1] - forceNodeVel[iTurb][iNodeForceLower][1] ); + cDriver_Output_to_FAST[iTurb].w[iNodeVel] = forceNodeVel[iTurb][iNodeForceLower][2] + rInterp * (forceNodeVel[iTurb][iNodeForceLower+1][2] - forceNodeVel[iTurb][iNodeForceLower][2] ); + } + } + } } void fast::OpenFAST::computeTorqueThrust(int iTurbGlob, std::vector & torque, std::vector & thrust) { @@ -684,7 +667,7 @@ void fast::OpenFAST::computeTorqueThrust(int iTurbGlob, std::vector & to for (int k=0; k < get_numBladesLoc(iTurbLoc); k++) { for (int j=0; j < numForcePtsBlade[iTurbLoc]; j++) { int iNode = 1 + numForcePtsBlade[iTurbLoc]*k + j ; - + thrust[0] = thrust[0] + cDriver_Input_from_FAST[iTurbLoc].fx[iNode] ; thrust[1] = thrust[1] + cDriver_Input_from_FAST[iTurbLoc].fy[iNode] ; thrust[2] = thrust[2] + cDriver_Input_from_FAST[iTurbLoc].fz[iNode] ; @@ -692,294 +675,272 @@ 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]; + + 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] + 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] ; - } } } fast::ActuatorNodeType fast::OpenFAST::getVelNodeType(int iTurbGlob, int iNode) { - // Return the type of velocity node for the given node number. The node ordering (from FAST) is - // Node 0 - Hub node - // Blade 1 nodes - // Blade 2 nodes - // Blade 3 nodes - // Tower nodes - - int iTurbLoc = get_localTurbNo(iTurbGlob); - for(int j=0; j < iTurbLoc; j++) iNode = iNode - get_numVelPtsLoc(iTurbGlob); - if (iNode) { - if ( (iNode + 1 - (get_numVelPts(iTurbLoc) - get_numVelPtsTwr(iTurbLoc)) ) > 0) { - return TOWER; - } - else { - return BLADE; + // Return the type of velocity node for the given node number. The node ordering (from FAST) is + // Node 0 - Hub node + // Blade 1 nodes + // Blade 2 nodes + // Blade 3 nodes + // Tower nodes + + int iTurbLoc = get_localTurbNo(iTurbGlob); + for(int j=0; j < iTurbLoc; j++) iNode = iNode - get_numVelPtsLoc(iTurbGlob); + if (iNode) { + if ( (iNode + 1 - (get_numVelPts(iTurbLoc) - get_numVelPtsTwr(iTurbLoc)) ) > 0 ) { + return TOWER; + } else { + return BLADE; + } + } else { + return HUB; } - } - else { - return HUB; - } - } fast::ActuatorNodeType fast::OpenFAST::getForceNodeType(int iTurbGlob, int iNode) { - // Return the type of actuator force node for the given node number. The node ordering (from FAST) is - // Node 0 - Hub node - // Blade 1 nodes - // Blade 2 nodes - // Blade 3 nodes - // Tower nodes - - int iTurbLoc = get_localTurbNo(iTurbGlob); - for(int j=0; j < iTurbLoc; j++) iNode = iNode - get_numForcePtsLoc(iTurbGlob); - if (iNode) { - if ( (iNode + 1 - (get_numForcePts(iTurbLoc) - get_numForcePtsTwr(iTurbLoc)) ) > 0) { - return TOWER; - } - else { - return BLADE; + // Return the type of actuator force node for the given node number. The node ordering (from FAST) is + // Node 0 - Hub node + // Blade 1 nodes + // Blade 2 nodes + // Blade 3 nodes + // Tower nodes + + int iTurbLoc = get_localTurbNo(iTurbGlob); + for(int j=0; j < iTurbLoc; j++) iNode = iNode - get_numForcePtsLoc(iTurbGlob); + if (iNode) { + if ( (iNode + 1 - (get_numForcePts(iTurbLoc) - get_numForcePtsTwr(iTurbLoc)) ) > 0 ) { + return TOWER; + } else { + return BLADE; + } + } else { + return HUB; } - } - else { - return HUB; - } - } void fast::OpenFAST::allocateMemory() { - - for (int iTurb=0; iTurb < nTurbinesGlob; iTurb++) { - if (dryRun) { - if(worldMPIRank == 0) { - std::cout << "iTurb = " << iTurb << " turbineMapGlobToProc[iTurb] = " << turbineMapGlobToProc[iTurb] << std::endl ; - } - } - if(worldMPIRank == turbineMapGlobToProc[iTurb]) { - turbineMapProcToGlob[nTurbinesProc] = iTurb; - reverseTurbineMapProcToGlob[iTurb] = nTurbinesProc; - nTurbinesProc++ ; + + for (int iTurb=0; iTurb < nTurbinesGlob; iTurb++) { + if (dryRun) { + if(worldMPIRank == 0) { + std::cout << "iTurb = " << iTurb << " turbineMapGlobToProc[iTurb] = " << turbineMapGlobToProc[iTurb] << std::endl ; + } + } + if(worldMPIRank == turbineMapGlobToProc[iTurb]) { + turbineMapProcToGlob[nTurbinesProc] = iTurb; + reverseTurbineMapProcToGlob[iTurb] = nTurbinesProc; + nTurbinesProc++ ; + } + turbineSetProcs.insert(turbineMapGlobToProc[iTurb]); } - turbineSetProcs.insert(turbineMapGlobToProc[iTurb]); - } - - int nProcsWithTurbines=0; - turbineProcs.resize(turbineSetProcs.size()); - - for (std::set::const_iterator p = turbineSetProcs.begin(); p != turbineSetProcs.end(); p++) { - turbineProcs[nProcsWithTurbines] = *p; - nProcsWithTurbines++ ; - } - - 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() ; + + int nProcsWithTurbines=0; + turbineProcs.resize(turbineSetProcs.size()); + + for (std::set::const_iterator p = turbineSetProcs.begin(); p != turbineSetProcs.end(); p++) { + turbineProcs[nProcsWithTurbines] = *p; + nProcsWithTurbines++ ; } - - } - - // Construct a group containing all procs running atleast 1 turbine in FAST - MPI_Group_incl(worldMPIGroup, nProcsWithTurbines, &turbineProcs[0], &fastMPIGroup) ; - int fastMPIcommTag = MPI_Comm_create(mpiComm, fastMPIGroup, &fastMPIComm); - if (MPI_COMM_NULL != fastMPIComm) { - MPI_Comm_rank(fastMPIComm, &fastMPIRank); - } - - TurbID.resize(nTurbinesProc); - TurbineBasePos.resize(nTurbinesProc); - FASTInputFileName.resize(nTurbinesProc); - CheckpointFileRoot.resize(nTurbinesProc); - nacelle_cd.resize(nTurbinesProc); - nacelle_area.resize(nTurbinesProc); - air_density.resize(nTurbinesProc); - numBlades.resize(nTurbinesProc); - numForcePtsBlade.resize(nTurbinesProc); - numForcePtsTwr.resize(nTurbinesProc); - numVelPtsBlade.resize(nTurbinesProc); - numVelPtsTwr.resize(nTurbinesProc); - forceNodeVel.resize(nTurbinesProc); - - for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { - - TurbineBasePos[iTurb].resize(3); - - int globProc = turbineMapProcToGlob[iTurb]; - TurbID[iTurb] = globTurbineData[globProc].TurbID; - FASTInputFileName[iTurb] = globTurbineData[globProc].FASTInputFileName ; - CheckpointFileRoot[iTurb] = globTurbineData[globProc].FASTRestartFileName ; - for(int i=0;i<3;i++) { - TurbineBasePos[iTurb][i] = globTurbineData[globProc].TurbineBasePos[i]; + 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() ; + } } - numForcePtsBlade[iTurb] = globTurbineData[globProc].numForcePtsBlade; - numForcePtsTwr[iTurb] = globTurbineData[globProc].numForcePtsTwr; - nacelle_cd[iTurb] = globTurbineData[globProc].nacelle_cd; - nacelle_area[iTurb] = globTurbineData[globProc].nacelle_area; - air_density[iTurb] = globTurbineData[globProc].air_density; - } + // Construct a group containing all procs running atleast 1 turbine in FAST + MPI_Group_incl(worldMPIGroup, nProcsWithTurbines, &turbineProcs[0], &fastMPIGroup) ; + int fastMPIcommTag = MPI_Comm_create(mpiComm, fastMPIGroup, &fastMPIComm); + if (MPI_COMM_NULL != fastMPIComm) { + MPI_Comm_rank(fastMPIComm, &fastMPIRank); + } - // Allocate memory for Turbine datastructure for all turbines - FAST_AllocateTurbines(&nTurbinesProc, &ErrStat, ErrMsg); - - // Allocate memory for OpFM Input types in FAST - cDriver_Input_from_FAST.resize(nTurbinesProc) ; - cDriver_Output_to_FAST.resize(nTurbinesProc) ; + TurbID.resize(nTurbinesProc); + TurbineBasePos.resize(nTurbinesProc); + FASTInputFileName.resize(nTurbinesProc); + CheckpointFileRoot.resize(nTurbinesProc); + nacelle_cd.resize(nTurbinesProc); + nacelle_area.resize(nTurbinesProc); + air_density.resize(nTurbinesProc); + numBlades.resize(nTurbinesProc); + numForcePtsBlade.resize(nTurbinesProc); + numForcePtsTwr.resize(nTurbinesProc); + numVelPtsBlade.resize(nTurbinesProc); + numVelPtsTwr.resize(nTurbinesProc); + forceNodeVel.resize(nTurbinesProc); - cDriverSC_Input_from_FAST.resize(nTurbinesProc) ; - cDriverSC_Output_to_FAST.resize(nTurbinesProc) ; + for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { + TurbineBasePos[iTurb].resize(3); + + int globProc = turbineMapProcToGlob[iTurb]; + TurbID[iTurb] = globTurbineData[globProc].TurbID; + FASTInputFileName[iTurb] = globTurbineData[globProc].FASTInputFileName ; + CheckpointFileRoot[iTurb] = globTurbineData[globProc].FASTRestartFileName ; + for(int i=0;i<3;i++) { + TurbineBasePos[iTurb][i] = globTurbineData[globProc].TurbineBasePos[i]; + } + numForcePtsBlade[iTurb] = globTurbineData[globProc].numForcePtsBlade; + numForcePtsTwr[iTurb] = globTurbineData[globProc].numForcePtsTwr; + nacelle_cd[iTurb] = globTurbineData[globProc].nacelle_cd; + nacelle_area[iTurb] = globTurbineData[globProc].nacelle_area; + air_density[iTurb] = globTurbineData[globProc].air_density; + } + + // Allocate memory for Turbine datastructure for all turbines + FAST_AllocateTurbines(&nTurbinesProc, &ErrStat, ErrMsg); + + // Allocate memory for OpFM Input types in FAST + cDriver_Input_from_FAST.resize(nTurbinesProc) ; + cDriver_Output_to_FAST.resize(nTurbinesProc) ; + + cDriverSC_Input_from_FAST.resize(nTurbinesProc) ; + cDriverSC_Output_to_FAST.resize(nTurbinesProc) ; } void fast::OpenFAST::allocateTurbinesToProcsSimple() { - - // Allocate turbines to each processor - round robin fashion - int nProcs ; - MPI_Comm_size(mpiComm, &nProcs); - for(int j = 0; j < nTurbinesGlob; j++) turbineMapGlobToProc[j] = j % nProcs ; - + // Allocate turbines to each processor - round robin fashion + int nProcs ; + MPI_Comm_size(mpiComm, &nProcs); + for(int j = 0; j < nTurbinesGlob; j++) turbineMapGlobToProc[j] = j % nProcs ; } 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++) { - FAST_End(&iTurb, &stopTheProgram); + // Deallocate types we allocated earlier + + if (nTurbinesProc > 0) closeVelocityDataFile(nt_global, velNodeDataFile); + + if ( !dryRun) { + bool stopTheProgram = false; + for (int iTurb=0; iTurb < nTurbinesProc; iTurb++) { + FAST_End(&iTurb, &stopTheProgram); + } + FAST_DeallocateTurbines(&ErrStat, ErrMsg); } - FAST_DeallocateTurbines(&ErrStat, ErrMsg); - } - - 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); + + 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::readVelocityData(int nTimesteps) { - - int nTurbines; - - 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); - } + int nTurbines; - // 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); + hid_t velDataFile = H5Fopen(("velDatafile." + std::to_string(worldMPIRank) + ".h5").c_str(), H5F_ACC_RDWR, H5P_DEFAULT); - 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] ); + { + hid_t attr = H5Aopen(velDataFile, "nTurbines", H5P_DEFAULT); + herr_t ret = H5Aread(attr, H5T_NATIVE_INT, &nTurbines) ; + H5Aclose(attr); } - herr_t status = H5Dclose(dset_id); + // 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); + } } 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); + 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); + } - { - 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 ; - 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 ; - 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); - 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); - 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); + } - 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); } - - } else { - // Open the file in append mode - velDataFile = H5Fopen(("velDatafile." + std::to_string(worldMPIRank) + ".h5").c_str(), H5F_ACC_RDWR, H5P_DEFAULT); - } - return velDataFile; + return velDataFile; } herr_t fast::OpenFAST::closeVelocityDataFile(int nt_global, hid_t velDataFile) { - - herr_t status = H5Fclose(velDataFile) ; - return status; + herr_t status = H5Fclose(velDataFile) ; + return status; } - void fast::OpenFAST::backupVelocityDataFile(int curTimeStep, hid_t & velDataFile) { closeVelocityDataFile(curTimeStep, velDataFile); @@ -996,152 +957,131 @@ void fast::OpenFAST::backupVelocityDataFile(int curTimeStep, hid_t & velDataFile 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()); + 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; - H5Dclose(dset_id); - H5Sclose(dspace_id); - H5Sclose(mspace_id); + std::vector tmpVelData; + tmpVelData.resize(nVelPts * 6); - hid_t attr_id = H5Aopen_by_name(h5File, ".", "nTimesteps", H5P_DEFAULT, H5P_DEFAULT); - herr_t status = H5Awrite(attr_id, H5T_NATIVE_INT, &iTimestep); - status = H5Aclose(attr_id); + 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()); -void fast::OpenFAST::applyVelocityData(int iPrestart, int iTurb, OpFM_OutputType_t cDriver_Output_to_FAST, std::vector & velData) { + H5Dclose(dset_id); + H5Sclose(dspace_id); + H5Sclose(mspace_id); - 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]; - } + hid_t attr_id = H5Aopen_by_name(h5File, ".", "nTimesteps", H5P_DEFAULT, H5P_DEFAULT); + herr_t status = H5Awrite(attr_id, H5T_NATIVE_INT, &iTimestep); + status = H5Aclose(attr_id); } -void fast::OpenFAST::loadSuperController(const fast::fastInputs & fi) { +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]; + } +} - if(fi.scStatus) { +void fast::OpenFAST::loadSuperController(const fast::fastInputs & fi) { - scStatus = fi.scStatus; - scLibFile = fi.scLibFile; + if(fi.scStatus) { - // open the library - scLibHandle = dlopen(scLibFile.c_str(), RTLD_LAZY); - if (!scLibHandle) { - std::cerr << "Cannot open library: " << dlerror() << '\n'; - } - - create_SuperController = (create_sc_t*) dlsym(scLibHandle, "create_sc"); - // reset errors - dlerror(); - const char *dlsym_error = dlerror(); - if (dlsym_error) { - std::cerr << "Cannot load symbol 'create_sc': " << dlsym_error << '\n'; - dlclose(scLibHandle); - } + scStatus = fi.scStatus; + scLibFile = fi.scLibFile; - destroy_SuperController = (destroy_sc_t*) dlsym(scLibHandle, "destroy_sc"); - // reset errors - dlerror(); - const char *dlsym_error_us = dlerror(); - if (dlsym_error_us) { - std::cerr << "Cannot load symbol 'destroy_sc': " << dlsym_error_us << '\n'; - dlclose(scLibHandle); - } + // open the library + scLibHandle = dlopen(scLibFile.c_str(), RTLD_LAZY); + if (!scLibHandle) { + std::cerr << "Cannot open library: " << dlerror() << '\n'; + } - sc = create_SuperController() ; + create_SuperController = (create_sc_t*) dlsym(scLibHandle, "create_sc"); + // reset errors + dlerror(); + const char *dlsym_error = dlerror(); + if (dlsym_error) { + std::cerr << "Cannot load symbol 'create_sc': " << dlsym_error << '\n'; + dlclose(scLibHandle); + } - numScInputs = fi.numScInputs; - numScOutputs = fi.numScOutputs; + destroy_SuperController = (destroy_sc_t*) dlsym(scLibHandle, "destroy_sc"); + // reset errors + dlerror(); + const char *dlsym_error_us = dlerror(); + if (dlsym_error_us) { + std::cerr << "Cannot load symbol 'destroy_sc': " << dlsym_error_us << '\n'; + dlclose(scLibHandle); + } - if ( (numScInputs > 0) && (numScOutputs > 0)) { - scOutputsGlob.resize(nTurbinesGlob*numScOutputs) ; - scInputsGlob.resize(nTurbinesGlob*numScInputs) ; - for (int iTurb=0; iTurb < nTurbinesGlob; iTurb++) { - for(int iInput=0; iInput < numScInputs; iInput++) { - scInputsGlob[iTurb*numScInputs + iInput] = 0.0 ; // Initialize to zero - } - for(int iOutput=0; iOutput < numScOutputs; iOutput++) { - scOutputsGlob[iTurb*numScOutputs + iOutput] = 0.0 ; // Initialize to zero - } - } + sc = create_SuperController() ; + + numScInputs = fi.numScInputs; + numScOutputs = fi.numScOutputs; + + if ( (numScInputs > 0) && (numScOutputs > 0)) { + scOutputsGlob.resize(nTurbinesGlob*numScOutputs) ; + scInputsGlob.resize(nTurbinesGlob*numScInputs) ; + for (int iTurb=0; iTurb < nTurbinesGlob; iTurb++) { + for(int iInput=0; iInput < numScInputs; iInput++) { + scInputsGlob[iTurb*numScInputs + iInput] = 0.0 ; // Initialize to zero + } + for(int iOutput=0; iOutput < numScOutputs; iOutput++) { + scOutputsGlob[iTurb*numScOutputs + iOutput] = 0.0 ; // Initialize to zero + } + } + } else { + std::cerr << "Make sure numScInputs and numScOutputs are greater than zero" << std::endl; + } } else { - std::cerr << "Make sure numScInputs and numScOutputs are greater than zero" << std::endl; + scStatus = false; + numScInputs = 0; + numScOutputs = 0; } - - } else { - scStatus = false; - numScInputs = 0; - numScOutputs = 0; - } - } - void fast::OpenFAST::fillScInputsGlob() { - - // Fills the global array containing inputs to the supercontroller from all turbines - for(int iTurb=0; iTurb < nTurbinesGlob; iTurb++) { - for(int iInput=0; iInput < numScInputs; iInput++) { - scInputsGlob[iTurb*numScInputs + iInput] = 0.0; // Initialize to zero + // Fills the global array containing inputs to the supercontroller from all turbines + + for(int iTurb=0; iTurb < nTurbinesGlob; iTurb++) { + for(int iInput=0; iInput < numScInputs; iInput++) { + scInputsGlob[iTurb*numScInputs + iInput] = 0.0; // Initialize to zero + } } - } - - for(int iTurb=0; iTurb < nTurbinesProc; iTurb++) { - for(int iInput=0; iInput < numScInputs; iInput++) { - scInputsGlob[turbineMapProcToGlob[iTurb]*numScInputs + iInput] = cDriverSC_Input_from_FAST[iTurb].toSC[iInput] ; + + for(int iTurb=0; iTurb < nTurbinesProc; iTurb++) { + for(int iInput=0; iInput < numScInputs; iInput++) { + scInputsGlob[turbineMapProcToGlob[iTurb]*numScInputs + iInput] = cDriverSC_Input_from_FAST[iTurb].toSC[iInput] ; + } } - } - - if (MPI_COMM_NULL != fastMPIComm) { - MPI_Allreduce(MPI_IN_PLACE, scInputsGlob.data(), numScInputs*nTurbinesGlob, MPI_DOUBLE, MPI_SUM, fastMPIComm) ; - } + if (MPI_COMM_NULL != fastMPIComm) { + MPI_Allreduce(MPI_IN_PLACE, scInputsGlob.data(), numScInputs*nTurbinesGlob, MPI_DOUBLE, MPI_SUM, fastMPIComm) ; + } } - void fast::OpenFAST::fillScOutputsLoc() { - - // Fills the local array containing outputs from the supercontroller to each turbine - - for(int iTurb=0; iTurb < nTurbinesProc; iTurb++) { - for(int iOutput=0; iOutput < numScOutputs; iOutput++) { - cDriverSC_Output_to_FAST[iTurb].fromSC[iOutput] = scOutputsGlob[turbineMapProcToGlob[iTurb]*numScOutputs + iOutput] ; + // Fills the local array containing outputs from the supercontroller to each turbine + for(int iTurb=0; iTurb < nTurbinesProc; iTurb++) { + for(int iOutput=0; iOutput < numScOutputs; iOutput++) { + cDriverSC_Output_to_FAST[iTurb].fromSC[iOutput] = scOutputsGlob[turbineMapProcToGlob[iTurb]*numScOutputs + iOutput] ; + } } - } - } - - - - - - - - - - - diff --git a/reg_tests/CMakeLists.txt b/reg_tests/CMakeLists.txt index 7f1beb912e..7101363587 100644 --- a/reg_tests/CMakeLists.txt +++ b/reg_tests/CMakeLists.txt @@ -37,6 +37,9 @@ option(CTEST_PLOT_ERRORS "Generate plots of regression test errors." OFF) # Set the OpenFAST executable configuration option and default set(CTEST_OPENFAST_EXECUTABLE "${CMAKE_BINARY_DIR}/glue-codes/openfast/openfast" CACHE FILEPATH "Specify the OpenFAST executable to use in testing.") +# Set the OpenFAST executable configuration option and default +set(CTEST_OPENFASTCPP_EXECUTABLE "${CMAKE_BINARY_DIR}/glue-codes/openfast-cpp/openfastcpp" CACHE FILEPATH "Specify the OpenFAST C++ executable to use in testing.") + # Set the AeroDyn executable configuration option and default set(CTEST_AERODYN_EXECUTABLE "${CMAKE_BINARY_DIR}/modules/aerodyn/aerodyn_driver" CACHE FILEPATH "Specify the AeroDyn driver executable to use in testing.") @@ -70,7 +73,7 @@ add_subdirectory("${CMAKE_CURRENT_LIST_DIR}/r-test") # build and seed the test directories with the data they need to run the tests file(MAKE_DIRECTORY ${CTEST_BINARY_DIR}) -foreach(regTest glue-codes/openfast modules/aerodyn modules/beamdyn modules/hydrodyn modules/subdyn) +foreach(regTest glue-codes/openfast glue-codes/openfast-cpp modules/aerodyn modules/beamdyn modules/hydrodyn modules/subdyn) file(MAKE_DIRECTORY ${CTEST_BINARY_DIR}/${regTest}) endforeach() diff --git a/reg_tests/CTestList.cmake b/reg_tests/CTestList.cmake index 23a8291278..ae727e08c2 100644 --- a/reg_tests/CTestList.cmake +++ b/reg_tests/CTestList.cmake @@ -84,6 +84,15 @@ function(of_regression_linear TESTNAME LABEL) regression(${TEST_SCRIPT} ${OPENFAST_EXECUTABLE} ${SOURCE_DIRECTORY} ${BUILD_DIRECTORY} ${TESTNAME} "${LABEL}") endfunction(of_regression_linear) +# openfast c-interface +function(of_regression_cpp TESTNAME LABEL) + set(TEST_SCRIPT "${CMAKE_CURRENT_LIST_DIR}/executeOpenfastCppRegressionCase.py") + set(OPENFAST_CPP_EXECUTABLE "${CTEST_OPENFASTCPP_EXECUTABLE}") + set(SOURCE_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}/..") + set(BUILD_DIRECTORY "${CTEST_BINARY_DIR}/glue-codes/openfast-cpp") + regression(${TEST_SCRIPT} ${OPENFAST_CPP_EXECUTABLE} ${SOURCE_DIRECTORY} ${BUILD_DIRECTORY} ${TESTNAME} "${LABEL}") +endfunction(of_regression_cpp) + # aerodyn function(ad_regression TESTNAME LABEL) set(TEST_SCRIPT "${CMAKE_CURRENT_LIST_DIR}/executeAerodynRegressionCase.py") @@ -155,6 +164,9 @@ of_regression("HelicalWake_OLAF" "openfast;aerodyn15;olaf" of_regression("EllipticalWing_OLAF" "openfast;aerodyn15;olaf") of_regression("StC_test_OC4Semi" "openfast;servodyn;hydrodyn;moordyn;offshore") +# OpenFAST C++ API test +of_regression_cpp("5MW_Land_DLL_WTurb_cpp" "openfast;cpp") + # AeroAcoustic regression test of_regression_aeroacoustic("IEA_LB_RWT-AeroAcoustics" "openfast;aerodyn15;aeroacoustics") diff --git a/reg_tests/executeOpenfastCppRegressionCase.py b/reg_tests/executeOpenfastCppRegressionCase.py new file mode 100644 index 0000000000..22a4bace4f --- /dev/null +++ b/reg_tests/executeOpenfastCppRegressionCase.py @@ -0,0 +1,174 @@ +# +# Copyright 2017 National Renewable Energy Laboratory +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +import os +import sys +basepath = os.path.sep.join(sys.argv[0].split(os.path.sep)[:-1]) if os.path.sep in sys.argv[0] else "." +sys.path.insert(0, os.path.sep.join([basepath, "lib"])) +import argparse +import shutil +import subprocess +import rtestlib as rtl +import openfastDrivers +import pass_fail +from errorPlotting import exportCaseSummary + +##### Helper functions +def ignoreBaselineItems(directory, contents): + itemFilter = ['linux-intel', 'linux-gnu', 'macos-gnu', 'windows-intel'] + caught = [] + for c in contents: + if c in itemFilter: + caught.append(c) + return tuple(caught) + +##### Main program + +### Store the python executable for future python calls +pythonCommand = sys.executable + +### Verify input arguments +parser = argparse.ArgumentParser(description="Executes OpenFAST and a regression test for a single test case.") +parser.add_argument("caseName", metavar="Case-Name", type=str, nargs=1, help="The name of the test case.") +parser.add_argument("executable", metavar="OpenFAST", type=str, nargs=1, help="The path to the OpenFAST executable.") +parser.add_argument("sourceDirectory", metavar="path/to/openfast_repo", type=str, nargs=1, help="The path to the OpenFAST repository.") +parser.add_argument("buildDirectory", metavar="path/to/openfast_repo/build", type=str, nargs=1, help="The path to the OpenFAST repository build directory.") +parser.add_argument("tolerance", metavar="Test-Tolerance", type=float, nargs=1, help="Tolerance defining pass or failure in the regression test.") +parser.add_argument("systemName", metavar="System-Name", type=str, nargs=1, help="The current system\'s name: [Darwin,Linux,Windows]") +parser.add_argument("compilerId", metavar="Compiler-Id", type=str, nargs=1, help="The compiler\'s id: [Intel,GNU]") +parser.add_argument("-p", "-plot", dest="plot", action='store_true', help="bool to include plots in failed cases") +parser.add_argument("-n", "-no-exec", dest="noExec", action='store_true', help="bool to prevent execution of the test cases") +parser.add_argument("-v", "-verbose", dest="verbose", action='store_true', help="bool to include verbose system output") + +args = parser.parse_args() + +caseName = args.caseName[0] +executable = os.path.abspath(args.executable[0]) +sourceDirectory = args.sourceDirectory[0] +buildDirectory = args.buildDirectory[0] +tolerance = args.tolerance[0] +systemName = args.systemName[0] +compilerId = args.compilerId[0] +plotError = args.plot +noExec = args.noExec +verbose = args.verbose + +# validate inputs +rtl.validateExeOrExit(executable) +rtl.validateDirOrExit(sourceDirectory) +if not os.path.isdir(buildDirectory): + os.makedirs(buildDirectory) + +### Map the system and compiler configurations to a solution set +# Internal names -> Human readable names +systemName_map = { + "darwin": "macos", + "linux": "linux", + "windows": "windows" +} +compilerId_map = { + "gnu": "gnu", + "intel": "intel" +} +# Build the target output directory name or choose the default +supportedBaselines = ["macos-gnu", "linux-intel", "linux-gnu", "windows-intel"] +targetSystem = systemName_map.get(systemName.lower(), "") +targetCompiler = compilerId_map.get(compilerId.lower(), "") +outputType = os.path.join(targetSystem+"-"+targetCompiler) +if outputType not in supportedBaselines: + outputType = supportedBaselines[0] +print("-- Using gold standard files with machine-compiler type {}".format(outputType)) + +### Build the filesystem navigation variables for running openfast on the test case +rtest = os.path.join(sourceDirectory, "reg_tests", "r-test") +moduleDirectory = os.path.join(rtest, "glue-codes", "openfast-cpp") +openfast_gluecode_directory = os.path.join(rtest, "glue-codes", "openfast") +inputsDirectory = os.path.join(moduleDirectory, caseName) +targetOutputDirectory = os.path.join(openfast_gluecode_directory, caseName.replace('_cpp', ''), outputType) +testBuildDirectory = os.path.join(buildDirectory, caseName) + +# verify all the required directories exist +if not os.path.isdir(rtest): + rtl.exitWithError("The test data directory, {}, does not exist. If you haven't already, run `git submodule update --init --recursive`".format(rtest)) +if not os.path.isdir(targetOutputDirectory): + rtl.exitWithError("The test data outputs directory, {}, does not exist. Try running `git submodule update`".format(targetOutputDirectory)) +if not os.path.isdir(inputsDirectory): + rtl.exitWithError("The test data inputs directory, {}, does not exist. Verify your local repository is up to date.".format(inputsDirectory)) + +# create the local output directory if it does not already exist +dst = os.path.join(buildDirectory, "5MW_Baseline") +src = os.path.join(openfast_gluecode_directory, "5MW_Baseline") +if not os.path.isdir(dst): + shutil.copytree(src, dst) +else: + names = os.listdir(src) + for name in names: + if name == "ServoData": + continue + srcname = os.path.join(src, name) + dstname = os.path.join(dst, name) + if os.path.isdir(srcname): + if not os.path.isdir(dstname): + shutil.copytree(srcname, dstname) + else: + shutil.copy2(srcname, dstname) + +if not os.path.isdir(testBuildDirectory): + shutil.copytree(inputsDirectory, testBuildDirectory, ignore=ignoreBaselineItems) + +### Run openfast on the test case +if not noExec: + cwd = os.getcwd() + os.chdir(testBuildDirectory) + print("** CWD: ", os.getcwd()) + caseInputFile = os.path.abspath("cDriver.yaml") + returnCode = openfastDrivers.runOpenfastCase(caseInputFile, executable) + if returnCode != 0: + rtl.exitWithError("") + os.chdir(cwd) + +### Build the filesystem navigation variables for running the regression test +localOutFile = os.path.join(testBuildDirectory, caseName + ".outb") +baselineOutFile = os.path.join(targetOutputDirectory, caseName.replace('_cpp', '') + ".outb") +rtl.validateFileOrExit(localOutFile) +rtl.validateFileOrExit(baselineOutFile) + +testData, testInfo, testPack = pass_fail.readFASTOut(localOutFile) +baselineData, baselineInfo, _ = pass_fail.readFASTOut(baselineOutFile) +performance = pass_fail.calculateNorms(testData, baselineData) +normalizedNorm = performance[:, 1] + +# export all case summaries +results = list(zip(testInfo["attribute_names"], [*performance])) +results_max = performance.max(axis=0) +exportCaseSummary(testBuildDirectory, caseName, results, results_max, tolerance) + +# failing case +if not pass_fail.passRegressionTest(normalizedNorm, tolerance): + if plotError: + from errorPlotting import finalizePlotDirectory, plotOpenfastError + for channel in testInfo["attribute_names"]: + try: + plotOpenfastError(localOutFile, baselineOutFile, channel) + except: + error = sys.exc_info()[1] + print("Error generating plots: {}".format(error)) + finalizePlotDirectory(localOutFile, testInfo["attribute_names"], caseName) + + sys.exit(1) + +# passing case +sys.exit(0) diff --git a/reg_tests/r-test b/reg_tests/r-test index ab933d3838..a48899cc79 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit ab933d383887ab4a0c369b0719fc2905fc859c79 +Subproject commit a48899cc79e16c97691d24d001790442db7f6f05