diff --git a/.github/workflows/CI_rosco-compile.yml b/.github/workflows/CI_rosco-compile.yml index 83cff4a2..97659a44 100644 --- a/.github/workflows/CI_rosco-compile.yml +++ b/.github/workflows/CI_rosco-compile.yml @@ -3,80 +3,146 @@ name: CI_rosco-compile # We run CI on push commits on all branches on: [push, pull_request] -# Specify FORTRAN compiler, used to be "gfortran-10" -env: - FORTRAN_COMPILER: gfortran # A workflow run is made up of one or more jobs that can run sequentially or in parallel jobs: - build: - name: Build (${{ matrix.os }}) - runs-on: ${{ matrix.os }} - strategy: - fail-fast: true - matrix: - os: ["ubuntu-latest", "windows-latest" , "macOS-latest"] - python-version: ["3.9"] - - steps: - - name: Checkout repository - uses: actions/checkout@v3 - - - name: Install conda/mamba - uses: conda-incubator/setup-miniconda@v2 - # https://github.com/marketplace/actions/setup-miniconda - with: - # To use mamba, uncomment here, comment out the miniforge line - mamba-version: "*" - # miniforge-version: "latest" - auto-update-conda: true - python-version: ${{ matrix.python-version }} - environment-file: environment.yml - activate-environment: test - auto-activate-base: false - miniforge-variant: Mambaforge - - # Install ROSCO toolbox - - name: Install ROSCO toolbox - shell: bash -l {0} - run: | - python setup.py install + build_pip: + name: Pip Build (${{ matrix.os }}) - ${{ matrix.python-version }} + runs-on: ${{ matrix.os }} + defaults: + run: + shell: bash -l {0} + + strategy: + fail-fast: false #true + matrix: + os: ["ubuntu-latest", "macOS-latest", "windows-latest"] + python-version: ["3.9", "3.10", "3.11"] + + steps: + - name: Setup GNU Fortran + # if: false == contains( matrix.os, 'windows') + uses: awvwgk/setup-fortran@v1 #modflowpy/install-intelfortran-action@v1 # + + - name: checkout repository + uses: actions/checkout@v4 + + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v4 + id: cp + with: + python-version: ${{ matrix.python-version }} + update-environment: true + + #- name: Setup tmate session + # if: contains( matrix.os, 'windows') + # uses: mxschmitt/action-tmate@v3 + + # Install ZeroMQ- should be done via Conda + #- name: Install ZeroMQ + # run: | + # sudo apt-get update + # sudo apt-get install libzmq3-dev + + - name: Pip Install ROSCO + run: | + '${{ steps.cp.outputs.python-path }}' -m pip install -e . - # Re-generate registry - - name: Generate Registry - shell: bash -l {0} - run: python rosco/controller/rosco_registry/write_registry.py + # Re-generate registry + - name: Generate Registry + run: | + '${{ steps.cp.outputs.python-path }}' rosco/controller/rosco_registry/write_registry.py + + - name: Test run + # skipping pip test on windows for now until wisdem pypi is ready + if: false == contains( matrix.os, 'windows') + run: | + cd Examples + '${{ steps.cp.outputs.python-path }}' 01_turbine_model.py + + + build_and_test_conda: + name: Conda Build (${{ matrix.os }}) - ${{ matrix.python-version }} + runs-on: ${{ matrix.os }} + defaults: + run: + shell: bash -el {0} + + strategy: + fail-fast: false #true + matrix: + os: ["ubuntu-latest", "macOS-latest", "windows-latest"] + python-version: ["3.9", "3.10", "3.11"] + + steps: + - name: checkout repository + uses: actions/checkout@v4 + + - uses: conda-incubator/setup-miniconda@v2 + # https://github.com/marketplace/actions/setup-miniconda + with: + #mamba-version: "*" + miniforge-version: "latest" + auto-update-conda: true + python-version: ${{ matrix.python-version }} + environment-file: environment.yml + activate-environment: test + auto-activate-base: false + + - name: Add dependencies windows specific + if: contains( matrix.os, 'windows') + run: | + conda install -y m2w64-toolchain libpython + gfortran --version + + - name: Add dependencies mac specific + if: contains( matrix.os, 'mac') + run: | + conda install -y compilers + gfortran --version + + # Install + - name: Debug + run: | + conda list + printenv + + #- name: Setup tmate session + # uses: mxschmitt/action-tmate@v3 + # with: + # detached: true + # if: contains( matrix.os, 'windows') + + - name: Conda Install ROSCO + run: | + python -m pip install -e . + + - name: Install OpenFAST + run: | + conda install openfast==3.5.1 - - name: Add dependencies windows - if: true == contains( matrix.os, 'windows') - run: | - conda install -y m2w64-toolchain - - - name: Add dependencies windows - if: true == contains( matrix.os, 'mac') - shell: bash -l {0} - run: | - conda install -y gfortran - - - name: Setup Workspace - run: | - cmake -E make_directory ${{runner.workspace}}/ROSCO/rosco/controller/build - - - name: Configure and Build - unix - if: false == contains( matrix.os, 'windows') - shell: bash -l {0} - working-directory: "${{runner.workspace}}/ROSCO/rosco/controller/build" - run: | - cmake \ - -DCMAKE_INSTALL_PREFIX:PATH=${{runner.workspace}}/ROSCO/rosco/controller/install \ - -DCMAKE_Fortran_COMPILER:STRING=${{env.FORTRAN_COMPILER}} \ - .. - cmake --build . --target install - - - name: Configure and Build - windows - if: true == contains( matrix.os, 'windows') - #shell: bash #-l {0} - working-directory: "${{runner.workspace}}/ROSCO/rosco/controller/build" - run: | - cmake -G "MinGW Makefiles" -DCMAKE_INSTALL_PREFIX="${{runner.workspace}}/ROSCO/rosco/controller/build" .. - cmake --build . --target install + - name: Generate Registry + run: | + python rosco/controller/rosco_registry/write_registry.py + + - name: Fix example file extensions + run: | + python Examples/Test_Cases/update_libdiscon_extension.py + + - name: Run ROSCO testing + run: | + cd rosco/test + python ROSCO_testing.py + + - name: Run regression testing + if: contains( matrix.os, 'ubuntu') + run: | + cd rosco/test + pytest . + + - name: Test walkthrough notebook + if: contains( matrix.os, 'ubuntu') + run: | + cd Examples + treon ROSCO_walkthrough.ipynb + diff --git a/.github/workflows/CI_rosco-pytools.yml b/.github/workflows/CI_rosco-pytools.yml deleted file mode 100644 index f2f3b6f5..00000000 --- a/.github/workflows/CI_rosco-pytools.yml +++ /dev/null @@ -1,226 +0,0 @@ -name: CI_rosco-pytools - -# We run CI on push commits on all branches -on: [push, pull_request] - -# Specify FORTRAN compiler, used to be "gfortran-10" -env: - FC: gfortran - -# A workflow run is made up of one or more jobs that can run sequentially or in parallel -jobs: - build: - name: Build (${{ matrix.os }}) - runs-on: ${{ matrix.os }} - strategy: - fail-fast: true - matrix: - os: ["ubuntu-latest", "windows-latest", "macOS-latest"] - python-version: ["3.9"] - defaults: - run: - shell: bash -l {0} - - steps: - - name: Checkout repository and submodules - uses: actions/checkout@v3 - with: - submodules: recursive - - - name: Install conda/mamba - uses: conda-incubator/setup-miniconda@v2 - # https://github.com/marketplace/actions/setup-miniconda - with: - # To use mamba, uncomment here, comment out the miniforge line - mamba-version: "*" - # miniforge-version: "latest" - auto-update-conda: true - python-version: ${{ matrix.python-version }} - environment-file: environment.yml - activate-environment: test - auto-activate-base: false - miniforge-variant: Mambaforge - - - # Install dependencies of ROSCO toolbox - - name: Add dependencies windows specific - if: true == contains( matrix.os, 'windows') - run: | - conda install -y m2w64-toolchain libpython - env: - FC: gfortran - - - - name: Add dependencies macOS specific - if: true == contains( matrix.os, 'macOS') - run: | - conda install compilers - - # Install ROSCO toolbox - - name: Install ROSCO toolbox on Windows - if: true == contains( matrix.os, 'windows') - env: - FC: gfortran - run: | - python setup.py develop --compile-rosco - - - name: Install ROSCO toolbox - if: false == contains( matrix.os, 'windows') - run: | - python setup.py develop --compile-rosco - - - run_examples: - name: Run Examples - runs-on: ${{ matrix.os }} - strategy: - fail-fast: true - matrix: - os: ["ubuntu-latest"] #, "macOS-latest"] - python-version: ["3.9"] - defaults: - run: - shell: bash -l {0} - - steps: - - name: Checkout repository and submodules - uses: actions/checkout@v3 - with: - submodules: recursive - - - name: Install conda/mamba - uses: conda-incubator/setup-miniconda@v2 - # https://github.com/marketplace/actions/setup-miniconda - with: - # To use mamba, uncomment here, comment out the miniforge line - mamba-version: "*" - # miniforge-version: "latest" - auto-update-conda: true - python-version: ${{ matrix.python-version }} - environment-file: environment.yml - activate-environment: test - auto-activate-base: false - miniforge-variant: Mambaforge - - # setup cmake - - name: Setup Workspace - run: cmake -E make_directory ${{runner.workspace}}/ROSCO/rosco/controller/build - - # Install dependencies of ROSCO toolbox - - name: Add dependencies ubuntu specific - run: | - conda install -y wisdem - - - name: Add pyFAST dependency - run: | - git clone http://github.com/OpenFAST/python-toolbox - cd python-toolbox - python -m pip install -e . - - # Install ZeroMQ - - name: Install ZeroMQ - run: | - sudo apt-get update - sudo apt-get install libzmq3-dev - - # Install ROSCO toolbox - - name: Install ROSCO toolbox - run: | - python setup.py develop - - - name: Configure and Build ROSCO - unix - working-directory: ${{runner.workspace}}/ROSCO/rosco/controller/build - run: | - cmake \ - -DCMAKE_INSTALL_PREFIX:PATH=${{runner.workspace}}/ROSCO/rosco/controller/install \ - -DZMQ_CLIENT=ON \ - -DCMAKE_Fortran_COMPILER:STRING=${{env.FORTRAN_COMPILER}} \ - .. - cmake --build . --target install - - # Install OpenFAST - - name: Install OpenFAST - run: | - conda install openfast==3.5 - - # Run examples - - name: Run examples - run: | - cd Examples - python test_examples.py - - # Test walkthrough notebook - - name: Test walkthrough notebook - run: | - cd Examples - treon ROSCO_walkthrough.ipynb - - # # Archive artifacts - # - name: Archive production artifacts - # if: success() || failure() - # uses: actions/upload-artifact@v3 - # with: - # name: ROSCO-artifacts - # path: | - # ${{runner.workspace}} - - run_testing: - name: Run Testing - runs-on: ${{ matrix.os }} - strategy: - fail-fast: true - matrix: - os: ["ubuntu-latest"] #, "macOS-latest"] - python-version: ["3.9"] - defaults: - run: - shell: bash -l {0} - - steps: - - name: Checkout repository and submodules - uses: actions/checkout@v3 - with: - submodules: recursive - - - name: Install conda/mamba - uses: conda-incubator/setup-miniconda@v2 - # https://github.com/marketplace/actions/setup-miniconda - with: - # To use mamba, uncomment here, comment out the miniforge line - mamba-version: "*" - # miniforge-version: "latest" - auto-update-conda: true - python-version: ${{ matrix.python-version }} - environment-file: environment.yml - activate-environment: test - auto-activate-base: false - miniforge-variant: Mambaforge - - - # Install dependencies of ROSCO toolbox - - name: Add dependencies - run: | - conda install -y wisdem - - - # Install ROSCO toolbox - - name: Install ROSCO toolbox - run: | - python setup.py install --compile-rosco - - # Install OpenFAST - - name: Install OpenFAST - run: | - conda install openfast==3.5 - - # Run ROSCO Testing - - name: Run ROSCO testing - run: | - cd rosco/test - python ROSCO_testing.py - - # Regression testing - - name: Run regression testing - run: | - cd ROSCO_testing - python test_checkpoint.py diff --git a/.github/workflows/Publish_ROSCO.yml b/.github/workflows/Publish_ROSCO.yml new file mode 100644 index 00000000..c6d7ba46 --- /dev/null +++ b/.github/workflows/Publish_ROSCO.yml @@ -0,0 +1,72 @@ +name: Build and upload to PyPI +# https://github.com/pypa/cibuildwheel/blob/main/examples/github-deploy.yml + +# Build on every branch push, tag push, and pull request change: +#on: [push, pull_request] +# Alternatively, to publish when a (published) GitHub Release is created, use the following: +on: + release: + types: + - published + +jobs: + build_wheels: + name: Build wheels on ${{ matrix.os }} + runs-on: ${{ matrix.os }} + strategy: + matrix: + os: [ubuntu-latest, windows-latest, macos-latest] + + steps: + - name: Setup GNU Fortran + uses: awvwgk/setup-fortran@v1 + + - name: Checkout + uses: actions/checkout@v4 + + - name: Build wheels + uses: pypa/cibuildwheel@v2.16.2 + + - uses: actions/upload-artifact@v4 + with: + name: cibw-wheels-${{ matrix.os }}-${{ strategy.job-index }} + path: ./wheelhouse/*.whl + + build_sdist: + name: Build source distribution + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v3 + + - name: Build sdist + run: pipx run build --sdist + + - uses: actions/upload-artifact@v4 + with: + name: cibw-sdist + path: dist/*.tar.gz + + upload_pypi: + needs: [build_wheels, build_sdist] + runs-on: ubuntu-latest + environment: pypi + permissions: + id-token: write # IMPORTANT: this permission is mandatory for trusted publishing + # upload to PyPI on every tag starting with 'v' + #if: github.event_name == 'push' && startsWith(github.ref, 'refs/tags/v') + # alternatively, to publish when a GitHub Release is created, use the following rule: + if: github.event_name == 'release' && github.event.action == 'published' + steps: + - uses: actions/download-artifact@v4 + with: + # unpacks all CIBW artifacts into dist/ + pattern: cibw-* + path: dist + merge-multiple: true + + - name: Upload to PyPI + uses: pypa/gh-action-pypi-publish@release/v1 + #with: + # user: __token__ + # password: ${{ secrets.pypi_password }} + # # To test: repository_url: https://test.pypi.org/legacy/ diff --git a/Examples/04_simple_sim.py b/Examples/04_simple_sim.py index b35d4312..6e1c469c 100644 --- a/Examples/04_simple_sim.py +++ b/Examples/04_simple_sim.py @@ -37,17 +37,19 @@ controller_params = inps['controller_params'] # Specify controller dynamic library path and name -this_dir = os.path.dirname(os.path.abspath(__file__)) -example_out_dir = os.path.join(this_dir,'examples_out') -if not os.path.isdir(example_out_dir): - os.makedirs(example_out_dir) + +#directories +rosco_dir = os.path.dirname(this_dir) +example_out_dir = os.path.join(this_dir,'examples_out') +os.makedirs(example_out_dir,exist_ok=True) if platform.system() == 'Windows': - lib_name = os.path.join(this_dir, '../rosco/controller/build/libdiscon.dll') + sfx = 'dll' elif platform.system() == 'Darwin': - lib_name = os.path.join(this_dir, '../rosco/controller/build/libdiscon.dylib') + sfx = 'dylib' else: - lib_name = os.path.join(this_dir, '../rosco/controller/build/libdiscon.so') + sfx = 'so' +lib_name = os.path.join(rosco_dir, 'lib', 'libdiscon.'+sfx) # # Load turbine model from saved pickle turbine = ROSCO_turbine.Turbine diff --git a/Examples/11_robust_tuning.py b/Examples/11_robust_tuning.py index 14e238d8..171968ae 100644 --- a/Examples/11_robust_tuning.py +++ b/Examples/11_robust_tuning.py @@ -3,7 +3,7 @@ Controller tuning to satisfy a robustness criteria ------------------------------------- NOTE: This example necessitates the mbc3 through either pyFAST or WEIS -pyFAST is the easiest to install by cloning https://github.com/OpenFAST/python-toolbox and +pyFAST is the easiest to install by cloning https://github.com/OpenFAST/openfast_toolbox and running `python setup.py develop` from your conda environment In this example: diff --git a/Examples/14_open_loop_control.py b/Examples/14_open_loop_control.py index 930e7f5d..608f86b6 100644 --- a/Examples/14_open_loop_control.py +++ b/Examples/14_open_loop_control.py @@ -90,13 +90,13 @@ ### Run OpenFAST using aeroelasticse tools -# Set rosco_dll if platform.system() == 'Windows': - rosco_dll = os.path.join(rosco_dir, 'rosco/controller/build/libdiscon.dll') + sfx = 'dll' elif platform.system() == 'Darwin': - rosco_dll = os.path.join(rosco_dir, 'rosco/controller/build/libdiscon.dylib') + sfx = 'dylib' else: - rosco_dll = os.path.join(rosco_dir, 'rosco/controller/build/libdiscon.so') + sfx = 'so' +rosco_dll = os.path.join(rosco_dir, 'lib', 'libdiscon.'+sfx) case_inputs = {} case_inputs[('ServoDyn','DLL_FileName')] = {'vals': [rosco_dll], 'group': 0} diff --git a/Examples/16_external_dll.py b/Examples/16_external_dll.py index 162694bb..909b1c5b 100644 --- a/Examples/16_external_dll.py +++ b/Examples/16_external_dll.py @@ -20,11 +20,12 @@ os.makedirs(example_out_dir,exist_ok=True) if platform.system() == 'Windows': - lib_name = os.path.realpath(os.path.join(this_dir, '../rosco/controller/build/libdiscon.dll')) + sfx = 'dll' elif platform.system() == 'Darwin': - lib_name = os.path.realpath(os.path.join(this_dir, '../rosco/controller/build/libdiscon.dylib')) + sfx = 'dylib' else: - lib_name = os.path.realpath(os.path.join(this_dir, '../rosco/controller/build/libdiscon.so')) + sfx = 'so' +lib_name = os.path.join(rosco_dir, 'lib', 'libdiscon.'+sfx) def main(): diff --git a/Examples/17a_zeromq_simple.py b/Examples/17a_zeromq_simple.py index 28b23cec..81a58541 100644 --- a/Examples/17a_zeromq_simple.py +++ b/Examples/17a_zeromq_simple.py @@ -22,12 +22,24 @@ import numpy as np import multiprocessing as mp -this_dir = os.path.dirname(os.path.abspath(__file__)) -example_out_dir = os.path.join(this_dir, "examples_out") TIME_CHECK = 30 DESIRED_YAW_OFFSET = 20 DESIRED_PITCH_OFFSET = np.deg2rad(2) * np.sin(0.1 * TIME_CHECK) + np.deg2rad(2) +#directories +this_dir = os.path.dirname(os.path.abspath(__file__)) +rosco_dir = os.path.dirname(this_dir) +example_out_dir = os.path.join(this_dir,'examples_out') +os.makedirs(example_out_dir,exist_ok=True) + +if platform.system() == 'Windows': + sfx = 'dll' +elif platform.system() == 'Darwin': + sfx = 'dylib' +else: + sfx = 'so' +lib_name = os.path.join(rosco_dir, 'lib', 'libdiscon.'+sfx) + def run_zmq(logfile=None): # Start the server at the following address network_address = "tcp://*:5555" @@ -75,19 +87,6 @@ def sim_rosco(): controller_params['ZMQ_Mode'] = 1 controller_params['ZMQ_UpdatePeriod'] = 0.025 - # Specify controller dynamic library path and name - this_dir = os.path.dirname(os.path.abspath(__file__)) - example_out_dir = os.path.join(this_dir, 'examples_out') - if not os.path.isdir(example_out_dir): - os.makedirs(example_out_dir) - - if platform.system() == 'Windows': - lib_name = os.path.join(this_dir, '../ROSCO/build/libdiscon.dll') - elif platform.system() == 'Darwin': - lib_name = os.path.join(this_dir, '../ROSCO/build/libdiscon.dylib') - else: - lib_name = os.path.join(this_dir, '../ROSCO/build/libdiscon.so') - # # Load turbine model from saved pickle turbine = ROSCO_turbine.Turbine turbine = turbine.load(os.path.join(example_out_dir, '01_NREL5MW_saved.p')) @@ -172,4 +171,4 @@ def sim_rosco(): p2 = mp.Process(target=sim_rosco) p2.start() p1.join() - p2.join() \ No newline at end of file + p2.join() diff --git a/Examples/18_pitch_offsets.py b/Examples/18_pitch_offsets.py index a30628f1..98700dea 100644 --- a/Examples/18_pitch_offsets.py +++ b/Examples/18_pitch_offsets.py @@ -21,11 +21,12 @@ os.makedirs(example_out_dir,exist_ok=True) if platform.system() == 'Windows': - lib_name = os.path.realpath(os.path.join(this_dir, '../ROSCO/build/libdiscon.dll')) + sfx = 'dll' elif platform.system() == 'Darwin': - lib_name = os.path.realpath(os.path.join(this_dir, '../ROSCO/build/libdiscon.dylib')) + sfx = 'dylib' else: - lib_name = os.path.realpath(os.path.join(this_dir, '../ROSCO/build/libdiscon.so')) + sfx = 'so' +lib_name = os.path.join(rosco_dir, 'lib', 'libdiscon.'+sfx) def main(): @@ -84,4 +85,4 @@ def main(): if __name__=="__main__": - main() \ No newline at end of file + main() diff --git a/Examples/20_active_wake_control.py b/Examples/20_active_wake_control.py index 795955aa..90b97d21 100644 --- a/Examples/20_active_wake_control.py +++ b/Examples/20_active_wake_control.py @@ -163,6 +163,7 @@ # Choose your implementation method AWC_Mode = 1 # 1 for SNL implementation, 2 for Coleman Transformation implementation + #directories this_dir = os.path.dirname(os.path.abspath(__file__)) rosco_dir = os.path.dirname(this_dir) @@ -170,11 +171,13 @@ os.makedirs(example_out_dir,exist_ok=True) if platform.system() == 'Windows': - lib_name = os.path.realpath(os.path.join(this_dir, '../rosco/controller/build/libdiscon.dll')) + sfx = 'dll' elif platform.system() == 'Darwin': - lib_name = os.path.realpath(os.path.join(this_dir, '../rosco/controller/build/libdiscon.dylib')) + sfx = 'dylib' else: - lib_name = os.path.realpath(os.path.join(this_dir, '../rosco/controller/build/libdiscon.so')) + sfx = 'so' +lib_name = os.path.join(rosco_dir, 'lib', 'libdiscon.'+sfx) + def main(): diff --git a/Examples/21_optional_inputs.py b/Examples/21_optional_inputs.py index 84e2efc6..3bb88895 100644 --- a/Examples/21_optional_inputs.py +++ b/Examples/21_optional_inputs.py @@ -11,22 +11,24 @@ import numpy as np +#directories this_dir = os.path.dirname(os.path.abspath(__file__)) rosco_dir = os.path.dirname(this_dir) - example_out_dir = os.path.join(this_dir,'examples_out') example_in_dir = os.path.join(this_dir,'example_inputs') +os.makedirs(example_out_dir,exist_ok=True) +if platform.system() == 'Windows': + sfx = 'dll' +elif platform.system() == 'Darwin': + sfx = 'dylib' +else: + sfx = 'so' +rosco_dll = os.path.join(rosco_dir, 'lib', 'libdiscon.'+sfx) -def main(): - # Set up rosco_dll - if platform.system() == 'Windows': - rosco_dll = os.path.join(rosco_dir, 'rosco/controller/build/libdiscon.dll') - elif platform.system() == 'Darwin': - rosco_dll = os.path.join(rosco_dir, 'rosco/controller/build/libdiscon.dylib') - else: - rosco_dll = os.path.join(rosco_dir, 'rosco/controller/build/libdiscon.so') + +def main(): # Load turbine model from saved pickle turbine = ROSCO_turbine.Turbine @@ -48,4 +50,4 @@ def main(): assert avi_fail == [0,-1], "Unexpected pass/fail" if __name__ == "__main__": - main() \ No newline at end of file + main() diff --git a/Examples/23_structural_control.py b/Examples/23_structural_control.py index 7e12a22d..7034c5cf 100644 --- a/Examples/23_structural_control.py +++ b/Examples/23_structural_control.py @@ -36,11 +36,13 @@ os.makedirs(example_out_dir,exist_ok=True) if platform.system() == 'Windows': - lib_name = os.path.realpath(os.path.join(this_dir, '../ROSCO/build/libdiscon.dll')) + sfx = 'dll' elif platform.system() == 'Darwin': - lib_name = os.path.realpath(os.path.join(this_dir, '../ROSCO/build/libdiscon.dylib')) + sfx = 'dylib' else: - lib_name = os.path.realpath(os.path.join(this_dir, '../ROSCO/build/libdiscon.so')) + sfx = 'so' +lib_name = os.path.join(rosco_dir, 'lib', 'libdiscon.'+sfx) + def main(): @@ -121,4 +123,4 @@ def main(): if __name__=="__main__": - main() \ No newline at end of file + main() diff --git a/Examples/24_floating_feedback.py b/Examples/24_floating_feedback.py index 6d802f93..00107001 100644 --- a/Examples/24_floating_feedback.py +++ b/Examples/24_floating_feedback.py @@ -34,11 +34,12 @@ os.makedirs(example_out_dir,exist_ok=True) if platform.system() == 'Windows': - lib_name = os.path.realpath(os.path.join(this_dir, '../rosco/controller/build/libdiscon.dll')) + sfx = 'dll' elif platform.system() == 'Darwin': - lib_name = os.path.realpath(os.path.join(this_dir, '../rosco/controller/build/libdiscon.dylib')) + sfx = 'dylib' else: - lib_name = os.path.realpath(os.path.join(this_dir, '../rosco/controller/build/libdiscon.so')) + sfx = 'so' +lib_name = os.path.join(rosco_dir, 'lib', 'libdiscon.'+sfx) def main(): @@ -162,4 +163,4 @@ def main(): if __name__=="__main__": - main() \ No newline at end of file + main() diff --git a/Examples/25_rotor_position_control.py b/Examples/25_rotor_position_control.py index fd0e0d2b..58b24f4f 100644 --- a/Examples/25_rotor_position_control.py +++ b/Examples/25_rotor_position_control.py @@ -16,18 +16,21 @@ import pandas as pd import matplotlib.pyplot as plt + #directories this_dir = os.path.dirname(os.path.abspath(__file__)) -rosco_dir = os.path.dirname(this_dir) +rosco_dir = os.path.dirname(this_dir) example_out_dir = os.path.join(this_dir,'examples_out') os.makedirs(example_out_dir,exist_ok=True) if platform.system() == 'Windows': - lib_name = os.path.realpath(os.path.join(this_dir, '../rosco/controller/build/libdiscon.dll')) + sfx = 'dll' elif platform.system() == 'Darwin': - lib_name = os.path.realpath(os.path.join(this_dir, '../rosco/controller/build/libdiscon.dylib')) + sfx = 'dylib' else: - lib_name = os.path.realpath(os.path.join(this_dir, '../rosco/controller/build/libdiscon.so')) + sfx = 'so' +lib_name = os.path.join(rosco_dir, 'lib', 'libdiscon.'+sfx) + def main(): @@ -117,4 +120,4 @@ def main(): if __name__=="__main__": - main() \ No newline at end of file + main() diff --git a/Examples/27_power_ref_control.py b/Examples/27_power_ref_control.py index 73831607..12d253da 100644 --- a/Examples/27_power_ref_control.py +++ b/Examples/27_power_ref_control.py @@ -31,11 +31,12 @@ os.makedirs(example_out_dir,exist_ok=True) if platform.system() == 'Windows': - lib_name = os.path.realpath(os.path.join(this_dir, '../rosco/controller/build/libdiscon.dll')) + sfx = 'dll' elif platform.system() == 'Darwin': - lib_name = os.path.realpath(os.path.join(this_dir, '../rosco/controller/build/libdiscon.dylib')) + sfx = 'dylib' else: - lib_name = os.path.realpath(os.path.join(this_dir, '../rosco/controller/build/libdiscon.so')) + sfx = 'so' +lib_name = os.path.join(rosco_dir, 'lib', 'libdiscon.'+sfx) def main(): @@ -105,4 +106,4 @@ def main(): r.run_FAST() if __name__=="__main__": - main() \ No newline at end of file + main() diff --git a/Examples/28_tower_resonance.py b/Examples/28_tower_resonance.py index b4fc3319..a06e66b0 100644 --- a/Examples/28_tower_resonance.py +++ b/Examples/28_tower_resonance.py @@ -26,11 +26,12 @@ os.makedirs(example_out_dir,exist_ok=True) if platform.system() == 'Windows': - lib_name = os.path.realpath(os.path.join(this_dir, '../ROSCO/build/libdiscon.dll')) + sfx = 'dll' elif platform.system() == 'Darwin': - lib_name = os.path.realpath(os.path.join(this_dir, '../ROSCO/build/libdiscon.dylib')) + sfx = 'dylib' else: - lib_name = os.path.realpath(os.path.join(this_dir, '../ROSCO/build/libdiscon.so')) + sfx = 'so' +libname = os.path.join(rosco_dir, 'lib', 'libdiscon.'+sfx) def main(): diff --git a/Examples/ROSCO_walkthrough.ipynb b/Examples/ROSCO_walkthrough.ipynb index 05afd0f6..9a262889 100644 --- a/Examples/ROSCO_walkthrough.ipynb +++ b/Examples/ROSCO_walkthrough.ipynb @@ -556,7 +556,7 @@ "Generator speed: 1173.7 RPM, Pitch angle: 11.6 deg, Power: 5000.0 kW, Est. wind Speed: 15.7 m/s\n", "Generator speed: 1173.7 RPM, Pitch angle: 11.6 deg, Power: 5000.0 kW, Est. wind Speed: 15.7 m/s\n", "Generator speed: 1173.7 RPM, Pitch angle: 11.6 deg, Power: 5000.0 kW, Est. wind Speed: 15.7 m/s\n", - "Shutting down ../rosco/controller/build/libdiscon.dylib\n" + "Shutting down ../lib/libdiscon.dylib\n" ] }, { @@ -580,7 +580,7 @@ "else:\n", " ext = 'so'\n", " \n", - "lib_name = (f'../rosco/controller/build/libdiscon.{ext}')\n", + "lib_name = (f'../lib/libdiscon.{ext}')\n", "\n", "# Load the simulator and controller interface\n", "controller_int = ROSCO_ci.ControllerInterface(lib_name,param_filename=param_file)\n", diff --git a/Examples/Test_Cases/5MW_Land_Simulink/NRELOffshrBsline5MW_Onshore_ServoDyn.dat b/Examples/Test_Cases/5MW_Land_Simulink/NRELOffshrBsline5MW_Onshore_ServoDyn.dat index 27bb021c..0cd26388 100644 --- a/Examples/Test_Cases/5MW_Land_Simulink/NRELOffshrBsline5MW_Onshore_ServoDyn.dat +++ b/Examples/Test_Cases/5MW_Land_Simulink/NRELOffshrBsline5MW_Onshore_ServoDyn.dat @@ -74,7 +74,7 @@ True GenTiStp - Method to stop the generator {T: timed using TimGen ---------------------- CABLE CONTROL ------------------------------------------- 0 CCmode - Cable control mode {0: none, 4: user-defined from Simulink/Labview, 5: user-defined from Bladed-style DLL} (switch) ---------------------- BLADED INTERFACE ---------------------------------------- [used only with Bladed Interface] -"../../ROSCO/build/libdiscon.dylib" DLL_FileName - Name/location of the dynamic library {.dll [Windows] or .so [Linux]} in the Bladed-DLL format (-) [used only with Bladed Interface] +"../../../lib/libdiscon.dylib" DLL_FileName - Name/location of the dynamic library {.dll [Windows] or .so [Linux]} in the Bladed-DLL format (-) [used only with Bladed Interface] "../NREL-5MW/DISCON.IN" DLL_InFile - Name of input file sent to the DLL (-) [used only with Bladed Interface] "DISCON" DLL_ProcName - Name of procedure in DLL to be called (-) [case sensitive; used only with DLL Interface] "default" DLL_DT - Communication interval for dynamic library (s) (or "default") [used only with Bladed Interface] diff --git a/Examples/Test_Cases/BAR_10/BAR_10_DISCON.IN b/Examples/Test_Cases/BAR_10/BAR_10_DISCON.IN index 8dc17b7d..ebbc56f5 100644 --- a/Examples/Test_Cases/BAR_10/BAR_10_DISCON.IN +++ b/Examples/Test_Cases/BAR_10/BAR_10_DISCON.IN @@ -1,9 +1,10 @@ ! Controller parameter input file for the BAR_10 wind turbine -! - File written using ROSCO version 2.8.0 controller tuning logic on 12/21/23 +! - File written using ROSCO version 2.8.0 controller tuning logic on 01/05/24 !------- SIMULATION CONTROL ------------------------------------------------------------ 1 ! LoggingLevel - {0: write no debug files, 1: write standard output .dbg-file, 2: LoggingLevel 1 + ROSCO LocalVars (.dbg2) 3: LoggingLevel 2 + complete avrSWAP-array (.dbg3)} 0 ! DT_Out - {Time step to output .dbg* files, or 0 to match sampling period of OpenFAST} +1 ! Ext_Interface - (0 - use standard bladed interface, 1 - Use the extened DLL interface introduced in OpenFAST 3.5.0.) 0 ! Echo - (0 - no Echo, 1 - Echo input data to .echo) !------- CONTROLLER FLAGS ------------------------------------------------- diff --git a/Examples/Test_Cases/BAR_10/BAR_10_ServoDyn.dat b/Examples/Test_Cases/BAR_10/BAR_10_ServoDyn.dat index b13a1212..9c06419a 100644 --- a/Examples/Test_Cases/BAR_10/BAR_10_ServoDyn.dat +++ b/Examples/Test_Cases/BAR_10/BAR_10_ServoDyn.dat @@ -74,7 +74,7 @@ True GenTiStp - Method to stop the generator {T: timed usin ---------------------- CABLE CONTROL ------------------------------------------- 0 CCmode - Cable control mode {0: none, 4: user-defined from Simulink/Labview, 5: user-defined from Bladed-style DLL} (switch) ---------------------- BLADED INTERFACE ---------------------------------------- [used only with Bladed Interface] -"../../ROSCO/install/lib/libdiscon.dylib" DLL_FileName - Name/location of the dynamic library {.dll [Windows] or .so [Linux]} in the Bladed-DLL format (-) [used only with Bladed Interface] +"../../../lib/libdiscon.dylib" DLL_FileName - Name/location of the dynamic library {.dll [Windows] or .so [Linux]} in the Bladed-DLL format (-) [used only with Bladed Interface] "BAR_10_DISCON.IN" DLL_InFile - Name of input file sent to the DLL (-) [used only with Bladed Interface] "DISCON" DLL_ProcName - Name of procedure in DLL to be called (-) [case sensitive; used only with DLL Interface] "default" DLL_DT - Communication interval for dynamic library (s) (or "default") [used only with Bladed Interface] diff --git a/Examples/Test_Cases/IEA-15-240-RWT-UMaineSemi/DISCON-UMaineSemi.IN b/Examples/Test_Cases/IEA-15-240-RWT-UMaineSemi/DISCON-UMaineSemi.IN index 1ec21c5c..cd58ccf8 100644 --- a/Examples/Test_Cases/IEA-15-240-RWT-UMaineSemi/DISCON-UMaineSemi.IN +++ b/Examples/Test_Cases/IEA-15-240-RWT-UMaineSemi/DISCON-UMaineSemi.IN @@ -1,9 +1,10 @@ ! Controller parameter input file for the IEA-15-240-RWT-UMaineSemi wind turbine -! - File written using ROSCO version 2.8.0 controller tuning logic on 12/21/23 +! - File written using ROSCO version 2.8.0 controller tuning logic on 01/05/24 !------- SIMULATION CONTROL ------------------------------------------------------------ 2 ! LoggingLevel - {0: write no debug files, 1: write standard output .dbg-file, 2: LoggingLevel 1 + ROSCO LocalVars (.dbg2) 3: LoggingLevel 2 + complete avrSWAP-array (.dbg3)} 0 ! DT_Out - {Time step to output .dbg* files, or 0 to match sampling period of OpenFAST} +1 ! Ext_Interface - (0 - use standard bladed interface, 1 - Use the extened DLL interface introduced in OpenFAST 3.5.0.) 0 ! Echo - (0 - no Echo, 1 - Echo input data to .echo) !------- CONTROLLER FLAGS ------------------------------------------------- diff --git a/Examples/Test_Cases/IEA-15-240-RWT-UMaineSemi/IEA-15-240-RWT-UMaineSemi_ServoDyn.dat b/Examples/Test_Cases/IEA-15-240-RWT-UMaineSemi/IEA-15-240-RWT-UMaineSemi_ServoDyn.dat index bb867d22..7ea5ee55 100644 --- a/Examples/Test_Cases/IEA-15-240-RWT-UMaineSemi/IEA-15-240-RWT-UMaineSemi_ServoDyn.dat +++ b/Examples/Test_Cases/IEA-15-240-RWT-UMaineSemi/IEA-15-240-RWT-UMaineSemi_ServoDyn.dat @@ -74,7 +74,7 @@ True GenTiStp - Method to stop the generator {T: timed usin ---------------------- CABLE CONTROL ------------------------------------------- 0 CCmode - Cable control mode {0: none, 4: user-defined from Simulink/Labview, 5: user-defined from Bladed-style DLL} (switch) ---------------------- BLADED INTERFACE ---------------------------------------- [used only with Bladed Interface] -"../../../rosco/controller/install/lib/libdiscon.so" DLL_FileName - Name/location of the dynamic library {.dll [Windows] or .so [Linux]} in the Bladed-DLL format (-) [used only with Bladed Interface] +"../../../lib/libdiscon.so" DLL_FileName - Name/location of the dynamic library {.dll [Windows] or .so [Linux]} in the Bladed-DLL format (-) [used only with Bladed Interface] "DISCON-UMaineSemi.IN" DLL_InFile - Name of input file sent to the DLL (-) [used only with Bladed Interface] "DISCON" DLL_ProcName - Name of procedure in DLL to be called (-) [case sensitive; used only with DLL Interface] "default" DLL_DT - Communication interval for dynamic library (s) (or "default") [used only with Bladed Interface] diff --git a/Examples/Test_Cases/MHK_RM1/MHK_RM1_DISCON.IN b/Examples/Test_Cases/MHK_RM1/MHK_RM1_DISCON.IN index 34d3d4a3..322213ff 100644 --- a/Examples/Test_Cases/MHK_RM1/MHK_RM1_DISCON.IN +++ b/Examples/Test_Cases/MHK_RM1/MHK_RM1_DISCON.IN @@ -1,9 +1,10 @@ ! Controller parameter input file for the MHK_RM1_Floating wind turbine -! - File written using ROSCO version 2.8.0 controller tuning logic on 12/21/23 +! - File written using ROSCO version 2.8.0 controller tuning logic on 01/05/24 !------- SIMULATION CONTROL ------------------------------------------------------------ 2 ! LoggingLevel - {0: write no debug files, 1: write standard output .dbg-file, 2: LoggingLevel 1 + ROSCO LocalVars (.dbg2) 3: LoggingLevel 2 + complete avrSWAP-array (.dbg3)} 0 ! DT_Out - {Time step to output .dbg* files, or 0 to match sampling period of OpenFAST} +1 ! Ext_Interface - (0 - use standard bladed interface, 1 - Use the extened DLL interface introduced in OpenFAST 3.5.0.) 0 ! Echo - (0 - no Echo, 1 - Echo input data to .echo) !------- CONTROLLER FLAGS ------------------------------------------------- diff --git a/Examples/Test_Cases/MHK_RM1/MHK_RM1_ServoDyn.dat b/Examples/Test_Cases/MHK_RM1/MHK_RM1_ServoDyn.dat index d7f43a12..c41f6a55 100644 --- a/Examples/Test_Cases/MHK_RM1/MHK_RM1_ServoDyn.dat +++ b/Examples/Test_Cases/MHK_RM1/MHK_RM1_ServoDyn.dat @@ -74,7 +74,7 @@ True GenTiStp - Method to stop the generator {T: timed using TimGen ---------------------- CABLE CONTROL ------------------------------------------- 0 CCmode - Cable control mode {0: none, 4: user-defined from Simulink/Labview, 5: user-defined from Bladed-style DLL} (switch) ---------------------- BLADED INTERFACE ---------------------------------------- [used only with Bladed Interface] -"../../../rosco/controller/install/lib/libdiscon.so" DLL_FileName - Name/location of the dynamic library {.dll [Windows] or .so [Linux]} in the Bladed-DLL format (-) [used only with Bladed Interface] +"../../../lib/libdiscon.so" DLL_FileName - Name/location of the dynamic library {.dll [Windows] or .so [Linux]} in the Bladed-DLL format (-) [used only with Bladed Interface] "MHK_RM1_DISCON.IN" DLL_InFile - Name of input file sent to the DLL (-) [used only with Bladed Interface] "DISCON" DLL_ProcName - Name of procedure in DLL to be called (-) [case sensitive; used only with DLL Interface] "default" DLL_DT - Communication interval for dynamic library (s) (or "default") [used only with Bladed Interface] diff --git a/Examples/Test_Cases/NREL-5MW/DISCON.IN b/Examples/Test_Cases/NREL-5MW/DISCON.IN index ee3b7b70..0d8024d1 100644 --- a/Examples/Test_Cases/NREL-5MW/DISCON.IN +++ b/Examples/Test_Cases/NREL-5MW/DISCON.IN @@ -1,9 +1,10 @@ ! Controller parameter input file for the NREL-5MW wind turbine -! - File written using ROSCO version 2.8.0 controller tuning logic on 12/21/23 +! - File written using ROSCO version 2.8.0 controller tuning logic on 01/05/24 !------- SIMULATION CONTROL ------------------------------------------------------------ 1 ! LoggingLevel - {0: write no debug files, 1: write standard output .dbg-file, 2: LoggingLevel 1 + ROSCO LocalVars (.dbg2) 3: LoggingLevel 2 + complete avrSWAP-array (.dbg3)} 0 ! DT_Out - {Time step to output .dbg* files, or 0 to match sampling period of OpenFAST} +1 ! Ext_Interface - (0 - use standard bladed interface, 1 - Use the extened DLL interface introduced in OpenFAST 3.5.0.) 0 ! Echo - (0 - no Echo, 1 - Echo input data to .echo) !------- CONTROLLER FLAGS ------------------------------------------------- diff --git a/Examples/Test_Cases/NREL-5MW/NRELOffshrBsline5MW_Onshore_ServoDyn.dat b/Examples/Test_Cases/NREL-5MW/NRELOffshrBsline5MW_Onshore_ServoDyn.dat index bc5e207a..41627843 100644 --- a/Examples/Test_Cases/NREL-5MW/NRELOffshrBsline5MW_Onshore_ServoDyn.dat +++ b/Examples/Test_Cases/NREL-5MW/NRELOffshrBsline5MW_Onshore_ServoDyn.dat @@ -74,7 +74,7 @@ True GenTiStp - Method to stop the generator {T: timed using TimGen ---------------------- CABLE CONTROL ------------------------------------------- 0 CCmode - Cable control mode {0: none, 4: user-defined from Simulink/Labview, 5: user-defined from Bladed-style DLL} (switch) ---------------------- BLADED INTERFACE ---------------------------------------- [used only with Bladed Interface] -"../../../rosco/controller/install/lib/libdiscon.so" DLL_FileName - Name/location of the dynamic library {.dll [Windows] or .so [Linux]} in the Bladed-DLL format (-) [used only with Bladed Interface] +"../../../lib/libdiscon.so" DLL_FileName - Name/location of the dynamic library {.dll [Windows] or .so [Linux]} in the Bladed-DLL format (-) [used only with Bladed Interface] "DISCON.IN" DLL_InFile - Name of input file sent to the DLL (-) [used only with Bladed Interface] "DISCON" DLL_ProcName - Name of procedure in DLL to be called (-) [case sensitive; used only with DLL Interface] "default" DLL_DT - Communication interval for dynamic library (s) (or "default") [used only with Bladed Interface] diff --git a/Examples/Test_Cases/NREL_2p8_127/NREL-2p8-127_DISCON.IN b/Examples/Test_Cases/NREL_2p8_127/NREL-2p8-127_DISCON.IN index acc3b083..31d0ed79 100644 --- a/Examples/Test_Cases/NREL_2p8_127/NREL-2p8-127_DISCON.IN +++ b/Examples/Test_Cases/NREL_2p8_127/NREL-2p8-127_DISCON.IN @@ -1,9 +1,10 @@ ! Controller parameter input file for the NREL-2p8-127 wind turbine -! - File written using ROSCO version 2.8.0 controller tuning logic on 12/21/23 +! - File written using ROSCO version 2.8.0 controller tuning logic on 01/05/24 !------- SIMULATION CONTROL ------------------------------------------------------------ 2 ! LoggingLevel - {0: write no debug files, 1: write standard output .dbg-file, 2: LoggingLevel 1 + ROSCO LocalVars (.dbg2) 3: LoggingLevel 2 + complete avrSWAP-array (.dbg3)} 0 ! DT_Out - {Time step to output .dbg* files, or 0 to match sampling period of OpenFAST} +1 ! Ext_Interface - (0 - use standard bladed interface, 1 - Use the extened DLL interface introduced in OpenFAST 3.5.0.) 0 ! Echo - (0 - no Echo, 1 - Echo input data to .echo) !------- CONTROLLER FLAGS ------------------------------------------------- diff --git a/Examples/Test_Cases/NREL_2p8_127/NREL-2p8-127_ServoDyn.dat b/Examples/Test_Cases/NREL_2p8_127/NREL-2p8-127_ServoDyn.dat index bdc2d74d..078bd87f 100644 --- a/Examples/Test_Cases/NREL_2p8_127/NREL-2p8-127_ServoDyn.dat +++ b/Examples/Test_Cases/NREL_2p8_127/NREL-2p8-127_ServoDyn.dat @@ -74,7 +74,7 @@ True GenTiStp - Method to stop the generator {T: timed usin ---------------------- CABLE CONTROL ------------------------------------------- 0 CCmode - Cable control mode {0: none, 4: user-defined from Simulink/Labview, 5: user-defined from Bladed-style DLL} (switch) ---------------------- BLADED INTERFACE ---------------------------------------- [used only with Bladed Interface] -"/pscratch/ndeveld/awc/ROSCO_B/ROSCO/build/libdiscon.so" DLL_FileName - Name/location of the dynamic library {.dll [Windows] or .so [Linux]} in the Bladed-DLL format (-) [used only with Bladed Interface] +"/pscratch/ndeveld/awc/ROSCO_B/lib/libdiscon.so" DLL_FileName - Name/location of the dynamic library {.dll [Windows] or .so [Linux]} in the Bladed-DLL format (-) [used only with Bladed Interface] "NREL-2p8-127_DISCON.IN" DLL_InFile - Name of input file sent to the DLL (-) [used only with Bladed Interface] "DISCON" DLL_ProcName - Name of procedure in DLL to be called (-) [case sensitive; used only with DLL Interface] default DLL_DT - Communication interval for dynamic library (s) (or "default") [used only with Bladed Interface] diff --git a/Examples/Test_Cases/update_libdiscon_extension.py b/Examples/Test_Cases/update_libdiscon_extension.py new file mode 100644 index 00000000..d85cb516 --- /dev/null +++ b/Examples/Test_Cases/update_libdiscon_extension.py @@ -0,0 +1,25 @@ +import glob +import platform +import os + +if platform.system() == 'Windows': + sfx = 'dll' +elif platform.system() == 'Darwin': + sfx = 'dylib' +else: + sfx = 'so' + +if __name__ == "__main__": + this_dir = os.path.dirname(os.path.abspath(__file__)) + servo_list = glob.glob(os.path.join(this_dir, '*/*Servo*.dat')) + + for ifile in servo_list: + # Read in current ServoDyn file + with open(ifile, "r") as f: + lines = f.readlines() + + # Write correction + with open(ifile, "w") as f: + for line in lines: + f.write(line.replace('libdiscon.so', f'libdiscon.{sfx}')) + diff --git a/Examples/Tune_Cases/IEA15MW_ExtInterface.yaml b/Examples/Tune_Cases/IEA15MW_ExtInterface.yaml index 0f65812d..32780386 100644 --- a/Examples/Tune_Cases/IEA15MW_ExtInterface.yaml +++ b/Examples/Tune_Cases/IEA15MW_ExtInterface.yaml @@ -54,6 +54,6 @@ controller_params: ps_percent: 0.8 # Percent peak shaving [%, <= 1 ], {default = 80%} DISCON: - DLL_FileName: /Users/dzalkind/Tools/ROSCO/ROSCO/build/libdiscon_copy.dylib + DLL_FileName: /Users/dzalkind/Tools/ROSCO/lib/libdiscon_copy.dylib DLL_InFile: /Users/dzalkind/Tools/ROSCO/Test_Cases/NREL-5MW/DISCON.IN - DLL_ProcName: DISCON \ No newline at end of file + DLL_ProcName: DISCON diff --git a/docs/source/api_change.rst b/docs/source/api_change.rst index 5bb95cdd..1e85bcad 100644 --- a/docs/source/api_change.rst +++ b/docs/source/api_change.rst @@ -9,18 +9,25 @@ The changes are tabulated according to the line number, and flag name. The line number corresponds to the resulting line number after all changes are implemented. Thus, be sure to implement each in order so that subsequent line numbers are correct. -2.8.0 to develop +2.8.0 to 2.9.0 ------------------------------- +**Flag to use exteneded Bladed Interface** + +* Set `Ext_Interface` to 1 to use the extened bladed interface with OpenFAST v3.5.0 and greater + **Gain scheduling of floating feedback** * The floating feedback gain can be scheduled on the low pass filtered wind speed signal. Note that Fl_Kp can now be an array. + +**Rotor position tracking** + +* Control the azimuth position of the rotor with `OL_Mode` of 2 using a PID torque controller with gains defined by `RP_Gains`. * Control all three blade pitch inputs in open loop -* Rotor position control of azimuth position with PID (RP_Gains) control inputs **New torque control mode settings** -* VS_ControlMode determines how the generator speed set point is determined: using the WSE or (P/K)^(1/3) -* VS_ConstPower determines whether constant power is used +* VS_ControlMode determines how the generator speed set point is determined: using the WSE (mode 2) or (P/K)^(1/3) (mode 3). The power signal in mode 3 is filtered using `VS_PwrFiltF`. +* VS_ConstPower determines whether constant power is used (0 is constant torque, 1 is constant power) **Multiple notch filters** @@ -56,31 +63,32 @@ New in ROSCO develop ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- Line Input Name Example Value ====== ======================= =============================================================================================================================================================================================================================================================== -13 VS_ConstPower 0 ! VS_ConstPower - Do constant power torque control, where above rated torque varies, 0 for constant torque} -17 PRC_Mode 0 ! PRC_Mode - Power reference tracking mode{0: use standard rotor speed set points, 1: use PRC rotor speed setpoints} -37 F_NumNotchFilts 1 ! F_NumNotchFilts - Number of notch filters placed on sensors -38 F_NotchFreqs 3.3550 ! F_NotchFreqs - Natural frequency of the notch filters. Array with length F_NumNotchFilts -39 F_NotchBetaNum 0.0000 ! F_NotchBetaNum - Damping value of numerator (determines the width of notch). Array with length F_NumNotchFilts, [-] -40 F_NotchBetaDen 0.2500 ! F_NotchBetaDen - Damping value of denominator (determines the depth of notch). Array with length F_NumNotchFilts, [-] -41 F_GenSpdNotch_N 0 ! F_GenSpdNotch_N - Number of notch filters on generator speed -42 F_GenSpdNotch_Ind 0 ! F_GenSpdNotch_Ind - Indices of notch filters on generator speed -43 F_TwrTopNotch_N 1 ! F_TwrTopNotch_N - Number of notch filters on tower top acceleration signal -44 F_TwrTopNotch_Ind 1 ! F_TwrTopNotch_Ind - Indices of notch filters on tower top acceleration signal -91 VS_PwrFiltF 0.3140 ! VS_PwrFiltF - Low pass filter on power used to determine generator speed set point. Only used in VS_ControlMode = 3. -97 PRC_Section !------- POWER REFERENCE TRACKING -------------------------------------- -98 PRC_n 2 ! PRC_n - Number of elements in PRC_WindSpeeds and PRC_GenSpeeds array -99 PRC_LPF_Freq 0.07854 ! PRC_LPF_Freq - Frequency of the low pass filter on the wind speed estimate used to set PRC_GenSpeeds [rad/s] -100 PRC_WindSpeeds 3.0000 25.0000 ! PRC_WindSpeeds - Array of wind speeds used in rotor speed vs. wind speed lookup table [m/s] -101 PRC_GenSpeeds 0.7917 0.7917 ! PRC_GenSpeeds - Array of generator speeds corresponding to PRC_WindSpeeds [rad/s] -102 Empty Line -127 TRA_ExclSpeed 0.00000 ! TRA_ExclSpeed - Rotor speed for exclusion [LSS, rad/s] -128 TRA_ExclBand 0.00000 ! TRA_ExclBand - Size of the rotor frequency exclusion band [LSS, rad/s]. Torque controller reference will be TRA_ExclSpeed +/- TRA_ExlBand/2 -129 TRA_RateLimit 0.00000e+00 ! TRA_RateLimit - Rate limit of change in rotor speed reference [LSS, rad/s]. Suggested to be VS_RefSpd/100. -144 Fl_n 1 ! Fl_n - Number of Fl_Kp gains in gain scheduling, optional with default of 1 -146 Fl_U 0.0000 ! Fl_U - Wind speeds for scheduling Fl_Kp, optional if Fl_Kp is single value [m/s] -160 Ind_Azimuth 0 ! Ind_Azimuth - The column in OL_Filename that contains the desired azimuth position in rad (used if OL_Mode = 2) -161 RP_Gains 0.0000 0.0000 0.0000 0.0000 ! RP_Gains - PID gains and Tf of derivative for rotor position control (used if OL_Mode = 2) -185 ZMQ_ID 0 ! ZMQ_ID - Integer identifier of turbine +7 Ext_Interface 1 ! Ext_Interface - (0 - use standard bladed interface, 1 - Use the extened DLL interface introduced in OpenFAST 3.5.0.) +14 VS_ConstPower 0 ! VS_ConstPower - Do constant power torque control, where above rated torque varies, 0 for constant torque} +18 PRC_Mode 0 ! PRC_Mode - Power reference tracking mode{0: use standard rotor speed set points, 1: use PRC rotor speed setpoints} +38 F_NumNotchFilts 1 ! F_NumNotchFilts - Number of notch filters placed on sensors +39 F_NotchFreqs 3.3550 ! F_NotchFreqs - Natural frequency of the notch filters. Array with length F_NumNotchFilts +40 F_NotchBetaNum 0.0000 ! F_NotchBetaNum - Damping value of numerator (determines the width of notch). Array with length F_NumNotchFilts, [-] +41 F_NotchBetaDen 0.2500 ! F_NotchBetaDen - Damping value of denominator (determines the depth of notch). Array with length F_NumNotchFilts, [-] +42 F_GenSpdNotch_N 0 ! F_GenSpdNotch_N - Number of notch filters on generator speed +43 F_GenSpdNotch_Ind 0 ! F_GenSpdNotch_Ind - Indices of notch filters on generator speed +44 F_TwrTopNotch_N 1 ! F_TwrTopNotch_N - Number of notch filters on tower top acceleration signal +45 F_TwrTopNotch_Ind 1 ! F_TwrTopNotch_Ind - Indices of notch filters on tower top acceleration signal +92 VS_PwrFiltF 0.3140 ! VS_PwrFiltF - Low pass filter on power used to determine generator speed set point. Only used in VS_ControlMode = 3. +98 PRC_Section !------- POWER REFERENCE TRACKING -------------------------------------- +99 PRC_n 2 ! PRC_n - Number of elements in PRC_WindSpeeds and PRC_GenSpeeds array +100 PRC_LPF_Freq 0.07854 ! PRC_LPF_Freq - Frequency of the low pass filter on the wind speed estimate used to set PRC_GenSpeeds [rad/s] +101 PRC_WindSpeeds 3.0000 25.0000 ! PRC_WindSpeeds - Array of wind speeds used in rotor speed vs. wind speed lookup table [m/s] +102 PRC_GenSpeeds 0.7917 0.7917 ! PRC_GenSpeeds - Array of generator speeds corresponding to PRC_WindSpeeds [rad/s] +103 Empty Line +128 TRA_ExclSpeed 0.00000 ! TRA_ExclSpeed - Rotor speed for exclusion [LSS, rad/s] +129 TRA_ExclBand 0.00000 ! TRA_ExclBand - Size of the rotor frequency exclusion band [LSS, rad/s]. Torque controller reference will be TRA_ExclSpeed +/- TRA_ExlBand/2 +130 TRA_RateLimit 0.00000e+00 ! TRA_RateLimit - Rate limit of change in rotor speed reference [LSS, rad/s]. Suggested to be VS_RefSpd/100. +145 Fl_n 1 ! Fl_n - Number of Fl_Kp gains in gain scheduling, optional with default of 1 +147 Fl_U 0.0000 ! Fl_U - Wind speeds for scheduling Fl_Kp, optional if Fl_Kp is single value [m/s] +161 Ind_Azimuth 0 ! Ind_Azimuth - The column in OL_Filename that contains the desired azimuth position in rad (used if OL_Mode = 2) +162 RP_Gains 0.0000 0.0000 0.0000 0.0000 ! RP_Gains - PID gains and Tf of derivative for rotor position control (used if OL_Mode = 2) +186 ZMQ_ID 0 ! ZMQ_ID - Integer identifier of turbine ====== ================= ====================================================================================================================================================================================================== ====== ================= ====================================================================================================================================================================================================== diff --git a/environment.yml b/environment.yml index d90d5165..fea8fbcd 100644 --- a/environment.yml +++ b/environment.yml @@ -3,16 +3,15 @@ channels: - defaults dependencies: + - cmake + - control - matplotlib - numpy - - pytest - - scipy - - pyYAML - pandas - - wisdem - - cmake - - pyparsing==2.4.7 - - control + - pyYAML + - pyparsing + - pytest - pyzmq - - mpi4py + - scipy - treon + - wisdem diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 00000000..ed0cf689 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,138 @@ +[build-system] +requires = ["setuptools", "cmake", "numpy", "wheel", "pyzmq"] +build-backend = "setuptools.build_meta" + +[project] +name = "rosco" +version = "2.9.0" +description = "A reference open source controller toolset for wind turbine applications." +readme = "README.md" +requires-python = ">=3.9" +license = {text = "Apache-2.0"} +keywords = ["wind", "turbine", "control", "", ""] +authors = [ + {name = "National Wind Technology Center, NREL", email = "daniel.zalkind@nrel.gov" } +] +maintainers = [ + {name = "Daniel Zalkind", email = "daniel.zalkind@nrel.gov" } +] +classifiers = [ # Optional + # How mature is this project? Common values are + # 3 - Alpha + # 4 - Beta + # 5 - Production/Stable + "Development Status :: 4 - Beta", + + # Indicate who your project is intended for + "Intended Audience :: Science/Research", + "Topic :: Scientific/Engineering", + + "License :: OSI Approved :: Apache Software License", + + # Specify the Python versions you support here. In particular, ensure + # that you indicate you support Python 3. These classifiers are *not* + # checked by "pip install". See instead "python_requires" below. + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.9", + "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", + "Programming Language :: Python :: 3 :: Only", + "Programming Language :: Fortran", +] + +dependencies = [ + "control", + "numpy", + "matplotlib", + "scipy", + "pandas", + "pyparsing", + "pyYAML", + "pyzmq", + "treon", + "wisdem", +] + +# List additional groups of dependencies here (e.g. development +# dependencies). Users will be able to install these using the "extras" +# syntax, for example: +# +# $ pip install sampleproject[dev] +# +# Similar to `dependencies` above, these must be valid existing +# projects. +[project.optional-dependencies] # Optional +dev = ["cmake"] +test = ["pytest"] + +# List URLs that are relevant to your project +# +# This field corresponds to the "Project-URL" and "Home-Page" metadata fields: +# https://packaging.python.org/specifications/core-metadata/#project-url-multiple-use +# https://packaging.python.org/specifications/core-metadata/#home-page-optional +# +# Examples listed include a pattern for specifying where the package tracks +# issues, where the source is hosted, where to say thanks to the package +# maintainers, and where to support the project financially. The key is +# what's used to render the link text on PyPI. +[project.urls] # Optional +"Homepage" = "https://github.com/NREL/ROSCO" +"Documentation" = "https://rosco.readthedocs.io" + +# This is configuration specific to the `setuptools` build backend. +# If you are using a different build backend, you will need to change this. +[tool.setuptools] +zip-safe = false +include-package-data = true + +#[tool.setuptools.packages] +#find = {} + +[tool.setuptools.packages.find] +#where = ["rosco"] +exclude = ["Matlab_Toolbox"] +namespaces = true + +[tool.setuptools.package-data] +# If there are data files included in your packages that need to be +# installed, specify them here. +"*" = ["*.yaml", "*.xlsx", "*.txt", "*.so", "*.lib", "*.pyd", "*.pdb", "*.dylib", "*.dll", "*.exe"] + +[tool.black] +line-length = 120 +target-version = ['py311'] +include = '\.pyi?$' +exclude = ''' +/( + \.git + | \.hg + | \.mypy_cache + | \.tox + | \.venv + | _build + | buck-out + | build + | dist +)/ +''' + +[tool.isort] +# https://github.com/PyCQA/isort +multi_line_output = "3" +include_trailing_comma = true +force_grid_wrap = false +use_parentheses = true +line_length = "120" +sections = ["FUTURE", "STDLIB", "THIRDPARTY", "FIRSTPARTY", "LOCALFOLDER"] +known_first_party = ["wisdem"] +length_sort = "1" +profile = "black" +skip_glob = ['__init__.py'] +atomic = true +#lines_after_imports = 2 +#lines_between_types = 1 +#src_paths=isort,test + +[tool.cibuildwheel] +skip = ["cp36-*", "cp37-*", "cp38-*", "*-win32"] +build-frontend = "build" diff --git a/rosco/controller/CMakeLists.txt b/rosco/controller/CMakeLists.txt index b450a30f..0d290851 100644 --- a/rosco/controller/CMakeLists.txt +++ b/rosco/controller/CMakeLists.txt @@ -1,5 +1,5 @@ cmake_minimum_required(VERSION 3.6) -project(ROSCO VERSION 2.8.0 LANGUAGES Fortran C) +project(ROSCO VERSION 2.9.0 LANGUAGES Fortran C) set(CMAKE_Fortran_MODULE_DIRECTORY "${CMAKE_BINARY_DIR}/ftnmods") @@ -26,13 +26,6 @@ set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -libs:static -free -static -fpp endif() endif() -# Enable ZMQ_Client if compiler flag is set -option(ZMQ_CLIENT "Enable use of ZeroMQ client" off) -if(ZMQ_Client) - set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -lzmq") -endif() - - set(SOURCES src/Constants.f90 src/ControllerBlocks.f90 @@ -69,11 +62,11 @@ else (NWTC_SYS_FILE) message(FATAL_ERROR "Cannot determine system file used with NWTC_Library") endif (NWTC_SYS_FILE) -# Library -if(ZMQ_CLIENT) +# ZMQ Library +find_package(PkgConfig) +pkg_check_modules(PC_ZeroMQ libzmq) +if(PC_ZeroMQ_FOUND) # Find ZeroMQ installation - find_package(PkgConfig) - pkg_check_modules(PC_ZeroMQ libzmq) find_path(ZeroMQ_INCLUDE_DIR NAMES zmq.h PATHS ${PC_ZeroMQ_INCLUDE_DIRS} @@ -84,20 +77,21 @@ if(ZMQ_CLIENT) ) include_directories(${ZeroMQ_INCLUDE_DIR}) + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -L${PC_ZeroMQ_LIBRARY_DIRS} -lzmq") # Compile C-based ZeroMQ client as object library add_compile_options(-I${ZeroMQ_INCLUDE_DIR} -l${ZeroMQ_LIBRARY} -fPIC) - add_library(zmq_client OBJECT ../src/zmq_client.c) + add_library(zmq_client OBJECT src/zmq_client.c) # Add definition add_definitions(-DZMQ_CLIENT="TRUE") - set(SOURCES ${SOURCES} + set(SOURCES ${SOURCES} $) endif() add_library(discon SHARED ${SOURCES}) -if(ZMQ_CLIENT) +if(PC_ZeroMQ_FOUND) target_include_directories(discon PUBLIC ${ZeroMQ_INCLUDE_DIR}) target_link_libraries(discon PUBLIC ${ZeroMQ_LIBRARY}) endif() @@ -114,5 +108,5 @@ set(linuxDefault ${CMAKE_INSTALL_PREFIX} STREQUAL "/usr/local") set(windowsDefault ${CMAKE_INSTALL_PREFIX} STREQUAL "C:\\Program Files (x86)") if(${linuxDefault} OR ${windowsDefault}) message("TRUE") - set(CMAKE_INSTALL_PREFIX "${CMAKE_SOURCE_DIR}/install") -endif() \ No newline at end of file + set(CMAKE_INSTALL_PREFIX "${CMAKE_SOURCE_DIR}/../../") +endif() diff --git a/rosco/controller/__init__.py b/rosco/controller/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/rosco/controller/rosco_registry/rosco_types.yaml b/rosco/controller/rosco_registry/rosco_types.yaml index a853772f..6f888f81 100644 --- a/rosco/controller/rosco_registry/rosco_types.yaml +++ b/rosco/controller/rosco_registry/rosco_types.yaml @@ -83,6 +83,9 @@ ControlParameters: Echo: <<: *integer description: 0 - no Echo, 1 - Echo input data to .echo + Ext_Interface: + <<: *integer + description: 0 - use standard bladed interface, 1 - Use the extened DLL interface introduced in OpenFAST 3.5.0. DT_Out: <<: *real diff --git a/rosco/controller/src/Constants.f90 b/rosco/controller/src/Constants.f90 index b85eefaa..56298879 100644 --- a/rosco/controller/src/Constants.f90 +++ b/rosco/controller/src/Constants.f90 @@ -14,7 +14,7 @@ MODULE Constants USE, INTRINSIC :: ISO_C_Binding - Character(*), PARAMETER :: rosco_version = 'v2.8.0' ! ROSCO version + Character(*), PARAMETER :: rosco_version = 'v2.9.0' ! ROSCO version INTEGER, PARAMETER :: DbKi = C_DOUBLE !< Default kind for double floating-point numbers INTEGER, PARAMETER :: ReKi = C_FLOAT !< Default kind for single floating-point numbers INTEGER, PARAMETER :: IntKi = C_INT !< Default kind for integer numbers diff --git a/rosco/controller/src/ROSCO_Types.f90 b/rosco/controller/src/ROSCO_Types.f90 index 6b8d00dd..48eaea0a 100644 --- a/rosco/controller/src/ROSCO_Types.f90 +++ b/rosco/controller/src/ROSCO_Types.f90 @@ -11,6 +11,7 @@ MODULE ROSCO_Types INTEGER(IntKi) :: ZMQ_ID ! 0000 - 9999, Identifier of the rosco, used for zeromq interface only INTEGER(IntKi) :: LoggingLevel ! 0 - write no debug files, 1 - write standard output .dbg-file, 2 - write standard output .dbg-file and complete avrSWAP-array .dbg2-file INTEGER(IntKi) :: Echo ! 0 - no Echo, 1 - Echo input data to .echo + INTEGER(IntKi) :: Ext_Interface ! 0 - use standard bladed interface, 1 - Use the extened DLL interface introduced in OpenFAST 3.5.0. REAL(DbKi) :: DT_Out ! Output time step INTEGER(IntKi) :: n_DT_Out ! output every this many steps INTEGER(IntKi) :: n_DT_ZMQ ! Send measurements to ZMQ after this many time steps diff --git a/rosco/controller/src/ReadSetParameters.f90 b/rosco/controller/src/ReadSetParameters.f90 index fe02e4db..fb159d2e 100644 --- a/rosco/controller/src/ReadSetParameters.f90 +++ b/rosco/controller/src/ReadSetParameters.f90 @@ -55,25 +55,29 @@ SUBROUTINE ReadAvrSWAP(avrSWAP, LocalVar, CntrPar) LocalVar%Azimuth = avrSWAP(60) LocalVar%NumBl = NINT(avrSWAP(61)) - ! Platform signals - LocalVar%PtfmTDX = avrSWAP(1001) - LocalVar%PtfmTDY = avrSWAP(1002) - LocalVar%PtfmTDZ = avrSWAP(1003) - LocalVar%PtfmRDX = avrSWAP(1004) - LocalVar%PtfmRDY = avrSWAP(1005) - LocalVar%PtfmRDZ = avrSWAP(1006) - LocalVar%PtfmTVX = avrSWAP(1007) - LocalVar%PtfmTVY = avrSWAP(1008) - LocalVar%PtfmTVZ = avrSWAP(1009) - LocalVar%PtfmRVX = avrSWAP(1010) - LocalVar%PtfmRVY = avrSWAP(1011) - LocalVar%PtfmRVZ = avrSWAP(1012) - LocalVar%PtfmTAX = avrSWAP(1013) - LocalVar%PtfmTAY = avrSWAP(1014) - LocalVar%PtfmTAZ = avrSWAP(1015) - LocalVar%PtfmRAX = avrSWAP(1016) - LocalVar%PtfmRAY = avrSWAP(1017) - LocalVar%PtfmRAZ = avrSWAP(1018) + if (CntrPar%Ext_Interface > 0) THEN ! Use extended bladed interface + + ! Platform signals + LocalVar%PtfmTDX = avrSWAP(1001) + LocalVar%PtfmTDY = avrSWAP(1002) + LocalVar%PtfmTDZ = avrSWAP(1003) + LocalVar%PtfmRDX = avrSWAP(1004) + LocalVar%PtfmRDY = avrSWAP(1005) + LocalVar%PtfmRDZ = avrSWAP(1006) + LocalVar%PtfmTVX = avrSWAP(1007) + LocalVar%PtfmTVY = avrSWAP(1008) + LocalVar%PtfmTVZ = avrSWAP(1009) + LocalVar%PtfmRVX = avrSWAP(1010) + LocalVar%PtfmRVY = avrSWAP(1011) + LocalVar%PtfmRVZ = avrSWAP(1012) + LocalVar%PtfmTAX = avrSWAP(1013) + LocalVar%PtfmTAY = avrSWAP(1014) + LocalVar%PtfmTAZ = avrSWAP(1015) + LocalVar%PtfmRAX = avrSWAP(1016) + LocalVar%PtfmRAY = avrSWAP(1017) + LocalVar%PtfmRAZ = avrSWAP(1018) + + ENDIF ! GenTemp = avrSWAP(1026) @@ -331,6 +335,7 @@ SUBROUTINE ReadControlParameterFileSub(CntrPar, LocalVar, accINFILE, accINFILE_s !----------------------- Simulation Control -------------------------- CALL ParseInput(FileLines,'LoggingLevel', CntrPar%LoggingLevel, accINFILE(1), ErrVar, .TRUE., UnEc=UnEc) CALL ParseInput(FileLines,'DT_Out', CntrPar%DT_Out, accINFILE(1), ErrVar, .TRUE., UnEc=UnEc) + CALL ParseInput(FileLines,'Ext_Interface', CntrPar%Ext_Interface, accINFILE(1), ErrVar, .TRUE., UnEc=UnEc) IF (ErrVar%aviFAIL < 0) RETURN !----------------- CONTROLLER FLAGS --------------------- @@ -362,7 +367,7 @@ SUBROUTINE ReadControlParameterFileSub(CntrPar, LocalVar, accINFILE, accINFILE_s !----------------- FILTER CONSTANTS --------------------- CALL ParseInput(FileLines, 'F_LPFCornerFreq', CntrPar%F_LPFCornerFreq, accINFILE(1), ErrVar, .FALSE., UnEc) CALL ParseInput(FileLines, 'F_LPFDamping', CntrPar%F_LPFDamping, accINFILE(1), ErrVar, CntrPar%F_LPFType == 1, UnEc) - CALL ParseInput(FileLines, 'F_NumNotchFilts', CntrPar%F_NumNotchFilts, accINFILE(1), ErrVar, .FALSE., UnEc) + CALL ParseInput(FileLines, 'F_NumNotchFilts', CntrPar%F_NumNotchFilts, accINFILE(1), ErrVar, .TRUE., UnEc) CALL ParseAry( FileLines, 'F_NotchFreqs', CntrPar%F_NotchFreqs, CntrPar%F_NumNotchFilts, accINFILE(1), ErrVar, CntrPar%F_NumNotchFilts == 0, UnEc) CALL ParseAry( FileLines, 'F_NotchBetaNum', CntrPar%F_NotchBetaNum, CntrPar%F_NumNotchFilts, accINFILE(1), ErrVar, CntrPar%F_NumNotchFilts == 0, UnEc) CALL ParseAry( FileLines, 'F_NotchBetaDen', CntrPar%F_NotchBetaDen, CntrPar%F_NumNotchFilts, accINFILE(1), ErrVar, CntrPar%F_NumNotchFilts == 0, UnEc) @@ -1399,6 +1404,13 @@ SUBROUTINE CheckInputs(LocalVar, CntrPar, avrSWAP, ErrVar, size_avcMSG) END IF IF (CntrPar%CC_Mode > 0) THEN + + ! Extended avrSWAP must be used + IF (CntrPar%Ext_Interface == 0) THEN + ErrVar%aviFAIL = -1 + ErrVar%ErrMsg = 'The OpenFAST extended bladed interface must be used with Ext_Interface > 0 in the DISCON' + ENDIF + IF (CntrPar%CC_ActTau .LE. 0) THEN ErrVar%aviFAIL = -1 ErrVar%ErrMsg = 'CC_ActTau must be greater than 0.' @@ -1413,6 +1425,14 @@ SUBROUTINE CheckInputs(LocalVar, CntrPar, avrSWAP, ErrVar, size_avcMSG) END IF IF (CntrPar%StC_Mode > 0) THEN + + ! Extended avrSWAP must be used + IF (CntrPar%Ext_Interface == 0) THEN + ErrVar%aviFAIL = -1 + ErrVar%ErrMsg = 'The OpenFAST extended bladed interface must be used with Ext_Interface > 0 in the DISCON' + ENDIF + + ! Check indices DO I = 1,CntrPar%StC_Group_N IF (CntrPar%StC_GroupIndex(I) < 2801) THEN ErrVar%aviFAIL = -1 diff --git a/rosco/test/ROSCO_testing.py b/rosco/test/ROSCO_testing.py index f938217b..f7967b03 100644 --- a/rosco/test/ROSCO_testing.py +++ b/rosco/test/ROSCO_testing.py @@ -13,17 +13,20 @@ ''' import numpy as np -import os, platform +import os +import platform import glob import multiprocessing as mp from rosco.toolbox.ofTools.fast_io.FAST_reader import InputReader_OpenFAST from rosco.toolbox.ofTools.case_gen.CaseGen_IEC import CaseGen_IEC from rosco.toolbox.ofTools.case_gen.runFAST_pywrapper import runFAST_pywrapper_batch -from matplotlib.backends.backend_pdf import FigureCanvasPdf, PdfPages +from matplotlib.backends.backend_pdf import PdfPages #, FigureCanvasPdf from rosco.toolbox.ofTools.fast_io import output_processing import matplotlib.pyplot as plt +this_dir = os.path.dirname(os.path.realpath(__file__)) +rosco_root = os.path.dirname( os.path.dirname( this_dir ) ) class ROSCO_testing(): @@ -35,7 +38,7 @@ def __init__(self, **kwargs): # Setup simulation parameters - self.runDir = os.path.join(os.path.dirname( os.path.realpath(__file__) ), 'testing' ) # directory to run simulations in + self.runDir = os.path.join(this_dir, 'testing' ) # directory to run simulations in self.wind_dir = None self.namebase = 'ROtest' # root name for output simulations self.FAST_exe = 'openfast_single' # name of openfast executable (may need full path) @@ -43,8 +46,8 @@ def __init__(self, **kwargs): self.FAST_ver = 'OpenFAST' # Fast version # Path to ROSCO controller - default to ROSCO Toolbox submodule try: - self.rosco_path = glob.glob(os.path.join(os.path.dirname(os.path.realpath(__file__)),'../rosco/controller/build/libdiscon.*'))[0] - except: + self.rosco_path = glob.glob(os.path.join(rosco_root, 'lib', 'libdiscon.*'))[0] + except Exception: print('No compiled ROSCO version found, please provide ROSCO_testing.rosco_path.') self.debug_level = 2 # debug level. 0 - no outputs, 1 - minimal outputs, 2 - all outputs self.overwrite = False # overwrite existing files? @@ -58,7 +61,7 @@ def __init__(self, **kwargs): # - Default to NREL 5MW self.Turbine_Class = 'I' self.Turbulence_Class = 'A' - self.FAST_directory = os.path.join(os.path.dirname(os.path.realpath(__file__)), '../../Examples/Test_Cases/NREL-5MW') + self.FAST_directory = os.path.join(rosco_root, 'Examples', 'Test_Cases','NREL-5MW') self.FAST_InputFile = 'NREL-5MW.fst' # Desired output channesl @@ -544,9 +547,6 @@ def print_results(self,outfiles): if __name__=='__main__': rt = ROSCO_testing() - this_dir = os.path.dirname(__file__) - - ## =================== INITIALIZATION =================== # Setup simulation parameters rt.namebase = 'IEA-15MW' # Base name for FAST files @@ -554,11 +554,13 @@ def print_results(self,outfiles): rt.Turbsim_exe = 'turbsim' # Turbsim executable path # path to compiled ROSCO controller if platform.system() == 'Windows': - rt.rosco_path = os.path.join(this_dir, '../rosco/controller/build/libdiscon.dll') + sfx = 'dll' elif platform.system() == 'Darwin': - rt.rosco_path = os.path.join(this_dir, '../rosco/controller/build/libdiscon.dylib') + sfx = 'dylib' else: - rt.rosco_path = os.path.join(this_dir, '../rosco/controller/build/libdiscon.so') + sfx = 'so' + rt.rosco_path = os.path.join(rosco_root, 'lib', 'libdiscon.'+sfx) + rt.debug_level = 2 # debug level. 0 - no outputs, 1 - minimal outputs, 2 - all outputs rt.overwrite = True # overwite fast sims? rt.cores = 1 # number of cores if multiprocessings @@ -569,7 +571,7 @@ def print_results(self,outfiles): # Setup turbine rt.Turbine_Class = 'I' rt.Turbulence_Class = 'B' - rt.FAST_directory = os.path.join(this_dir, '../../Examples/Test_Cases/IEA-15-240-RWT-UMaineSemi') + rt.FAST_directory = os.path.join(rosco_root, 'Examples','Test_Cases','IEA-15-240-RWT-UMaineSemi') rt.FAST_InputFile = 'IEA-15-240-RWT-UMaineSemi.fst' # Additional inputs diff --git a/rosco/test/run_Testing.py b/rosco/test/run_Testing.py index f10c0d01..314759cc 100644 --- a/rosco/test/run_Testing.py +++ b/rosco/test/run_Testing.py @@ -106,7 +106,7 @@ def run_testing(turbine2test, testtype, rosco_binaries=[], discon_files=[], **kw testtype = 'heavy' # lite, heavy, binary-comp, discon-comp # Only fill one of these if comparing controllers - rosco_binaries = [glob.glob(os.path.join(this_dir,'../ROSCO/build/libdiscon.*'))[0]] # Differently named libdiscons to compare + rosco_binaries = [glob.glob(os.path.join(this_dir,'..','..','lib','libdiscon.*'))[0]] # Differently named libdiscons to compare discon_files = [] # Differently named DISCON.IN files to compare # Run testing diff --git a/rosco/test/test_checkpoint.py b/rosco/test/test_checkpoint.py index 58abf1b0..c7d27a80 100644 --- a/rosco/test/test_checkpoint.py +++ b/rosco/test/test_checkpoint.py @@ -27,11 +27,12 @@ class RegressionTesting(unittest.TestCase): def test_restart(self): this_dir = os.path.dirname(os.path.abspath(__file__)) - rosco_dir = os.path.dirname(this_dir) + rosco_dir = os.path.dirname( os.path.dirname(this_dir) ) test_out_dir = os.path.join(this_dir, 'test_out') # Load yaml file (Open Loop Case) - parameter_filename = os.path.join(rosco_dir, 'Tune_Cases/IEA15MW.yaml') + tune_directory = os.path.join(rosco_dir,'Examples','Tune_Cases') + parameter_filename = os.path.join(tune_directory, 'IEA15MW.yaml') inps = load_rosco_yaml(parameter_filename) path_params = inps['path_params'] @@ -40,11 +41,12 @@ def test_restart(self): # Set rosco_dll if platform.system() == 'Windows': - rosco_dll = os.path.join(rosco_dir, 'ROSCO/build/libdiscon.dll') + sfx = 'dll' elif platform.system() == 'Darwin': - rosco_dll = os.path.join(rosco_dir, 'ROSCO/build/libdiscon.dylib') + sfx = 'dylib' else: - rosco_dll = os.path.join(rosco_dir, 'ROSCO/build/libdiscon.so') + sfx = 'so' + rosco_dll = os.path.join(rosco_dir, 'lib', 'libdiscon.'+sfx) case_inputs = {} case_inputs[('Fst', 'TMax')] = {'vals': [3.], 'group': 0} @@ -67,7 +69,8 @@ def test_restart(self): reader = InputReader_OpenFAST() writer = InputWriter_OpenFAST() reader.FAST_InputFile = path_params['FAST_InputFile'] - reader.FAST_directory = os.path.realpath(os.path.join( rosco_dir, 'Tune_Cases', path_params['FAST_directory'])) + reader.FAST_directory = os.path.realpath(os.path.join( tune_directory, path_params['FAST_directory'])) + reader.execute() writer.fst_vt = reader.fst_vt writer.FAST_runDirectory = test_out_dir @@ -115,6 +118,7 @@ def test_restart(self): fig, ax = op.plot_fast_out(cases=cases, showplot=False) plt.show() + print(fastout[1]['GenPwr'], fastout[0]['GenPwr']) self.check_relative_error(fastout[1]['GenPwr'], fastout[0]['GenPwr'], 1e-3) def check_relative_error(self, meas, real, tol): diff --git a/Examples/test_examples.py b/rosco/test/test_examples.py similarity index 89% rename from Examples/test_examples.py rename to rosco/test/test_examples.py index c182215e..2aa44cc0 100644 --- a/Examples/test_examples.py +++ b/rosco/test/test_examples.py @@ -36,7 +36,9 @@ ] def execute_script(fscript): - examples_dir = os.path.dirname(os.path.realpath(__file__)) + this_dir = os.path.dirname(os.path.abspath(__file__)) + rosco_dir = os.path.dirname( os.path.dirname(this_dir) ) + examples_dir = os.path.join(rosco_dir,'Examples') test_case_dir = os.path.join(examples_dir,'Test_Cases') # Go to location due to relative path use for airfoil files diff --git a/rosco/toolbox/__init__.py b/rosco/toolbox/__init__.py index 5b8f336b..3aa708bc 100644 --- a/rosco/toolbox/__init__.py +++ b/rosco/toolbox/__init__.py @@ -3,4 +3,4 @@ __author__ = """Nikhar J. Abbas and Daniel S. Zalkind""" __email__ = 'daniel.zalkind@nrel.gov' -__version__ = '2.8.0' +__version__ = '2.9.0' diff --git a/rosco/toolbox/inputs/toolbox_schema.yaml b/rosco/toolbox/inputs/toolbox_schema.yaml index 2d90620a..44573ae1 100644 --- a/rosco/toolbox/inputs/toolbox_schema.yaml +++ b/rosco/toolbox/inputs/toolbox_schema.yaml @@ -558,6 +558,12 @@ properties: type: number description: Time step to output .dbg* files, or 0 to match sampling period of OpenFAST default: 0 + Ext_Interface: + type: number + description: 0 - use standard bladed interface, 1 - Use the extened DLL interface introduced in OpenFAST 3.5.0. + minimum: 0 + maximum: 1 + default: 1 F_LPFType: type: number description: 1- first-order low-pass filter, 2- second-order low-pass filter (currently filters generator speed and pitch control signals diff --git a/rosco/toolbox/linear/getMats.py b/rosco/toolbox/linear/getMats.py new file mode 100644 index 00000000..7682b368 --- /dev/null +++ b/rosco/toolbox/linear/getMats.py @@ -0,0 +1,603 @@ +############################################################# +# Extract OpenFAST Matrices from linearization files # +# Authors: Srinivasa B. Ramisetti # +# Created: 01-July-2020 # +# E-mail: ramisettisrinivas@yahoo.com # +# Web: http://ramisetti.github.io # +############################################################# +#!/usr/bin/env python + +import os +from itertools import islice +import numpy as np +import re + +def _isbool(str): + flag=0 + if str=='T' or str=='True': + flag=1 + return flag + + +def findBladeTriplets(rotFrame, Desc, verbose=True): + + # Find the number of, and indices for, triplets in the rotating frame: + chkStr = [r'[Bb]lade \d', r'[Bb]lade [Rr]oot \d', r'BD_\d', r'[Bb]\d', r'[Bb]lade\d', r'PitchBearing\d', r'\d'] + + origDesc = Desc; + + # hack for ElastoDyn state names (remove unnecessary text in parenthesis) + for i in range(len(rotFrame)): + if Desc[i] is None: + raise Exception('Description not defined, likely a bug.') + ix = Desc[i].find('(internal DOF index = ') + if ix>0: + ix2 = Desc[i].find(')') + Desc[i] = Desc[i][:ix]+Desc[i][ix2+1:] + + + NTriplets = 0; # first initialize to zero + Triplets = []; + for i in range(len(rotFrame)): # loop through inputs/outputs/states + if rotFrame[i] == 'T': # this is in the rotating frame + Tmp = -1*np.ones(3) + foundTriplet = False + foundBladeNumber = False + for chk in chkStr: + BldNoCol = re.search(chk,Desc[i]) + if BldNoCol!=None: + foundBladeNumber = True + + Bldstr=BldNoCol.group() + # create another regular expression to find the + # exact match on a different blade: + strng = re.split(Bldstr,Desc[i],1) #this should return the strings before and after the match + + + FirstStr = strng[0] + Bldstr[:len(Bldstr)-1] + '.' + checkThisStr = FirstStr + strng[1] + + #we need to get rid of the special characters that + #may exist in Desc{}: + checkThisStr=checkThisStr.replace(')',r'\)').replace('(', r'\(').replace('^',r'\^') + FirstStr = FirstStr. replace(')',r'\)').replace('(', r'\(').replace('^',r'\^') + + k = int(Bldstr[len(Bldstr)-1]) + Tmp[k-1] = int(i) + break + + #print(Tmp,j) + + # find the other match values + if foundBladeNumber: + for j in range((i+1),len(rotFrame)): # loop through all remaining control inputs + #print(i, j) + if rotFrame[j]: # this is in the rotating frame + BldNoCol = re.search(checkThisStr, Desc[j]) + if BldNoCol!=None: + Num = re.search(FirstStr,Desc[j]).group(); + k = int(Num[len(Num)-1]); + Tmp[k-1] = int(j); # save the indices for the remaining blades + #TmpTmp=Tmp+1 + #print(BldNoCol.group(),i,j,k) + if ( (Tmp>-1).all() ): # true if all the elements of Tmp are nonzero; thus, we found a triplet of rotating indices + foundTriplet = True; + + Triplets.append(Tmp); # these are the indices for control input triplets in the rotating frame + + NTriplets = NTriplets + 1; # this is the number of control input triplets in the rotating frame + + # we'll set rotFrame to false so that we don't have to check the found channels again; also allows us to throw error if we have a rotating channel that doesn't have a unique match + for idx in Tmp: + id=int(idx) + rotFrame[id] = 0; + + break; + + if foundTriplet==False: + if verbose: + print('Rotating channel "', i, Desc[i], '" does not form a unique blade triplet. Blade(s) not found: ', np.array(np.where(Tmp == -1))+1 ) + else: + if verbose: + print( 'Could not find blade number in rotating channel "', Desc[i], '".') + #print(NTriplets) + #print(Triplets) + + return Triplets, NTriplets + +def reOrderByIdx_1D(arry,Indi): + tmp=[None]*len(arry) + j=0 + for i in Indi: + tmp[i]=arry[j] + j=j+1 + return tmp + +def reOrderByIdx_2D(arry,Indi,Indj): + #tmp=[[None]*arry.shape[0]]*arry.shape[1] + tmp=np.empty((arry.shape[0], arry.shape[1])) + #print(arry.shape, len(Indi), len(Indj)) + if arry.shape[0]!= len(Indi): + Indi=np.arange(0,arry.shape[0]) + if arry.shape[1]!= len(Indj): + Indj=np.arange(0,arry.shape[1]) + + p=0 + q=0 + for i in Indi: + q=0 + for j in Indj: + tmp[i][j]=arry[p][q] + q=q+1 + p=p+1 + return tmp + +def getStateOrderingIndx(matData): + + #StateOrderingIndx={} + StateOrderingIndx = np.arange(0,matData['NumStates']) + lastModName = ''; + lastModOrd = 0; + mod_nDOFs = 0; # number of DOFs in each module + sum_nDOFs2 = 0; # running total of second-order DOFs + sum_nDOFs1 = 0; # running total of first-order DOFs + indx_start = 0; # starting index of the modules + + for i in range(0,matData['NumStates']): + + tmp=(matData['DescStates'][i]); # name of the module whose states we are looking at + modName = tmp.split(' ')[0] + ModOrd = matData['StateDerivOrder'][i] + + if lastModName!=modName or lastModOrd!= ModOrd: + # this is the start of a new set of DOFs, so we'll set the + # indices for the last matrix + if lastModOrd == 2: + mod_nDOFs = int(mod_nDOFs/2); + #print('mmmm ', mod_nDOFs, indx_start, matData['ndof2'], i) + StateOrderingIndx[ indx_start :(indx_start+mod_nDOFs)] = sum_nDOFs2 + np.arange(0,mod_nDOFs); # q2 starts at 1 + StateOrderingIndx[ (indx_start+mod_nDOFs):(i)] = sum_nDOFs2 + matData['ndof2'] + np.arange(0,mod_nDOFs); # q2_dot starts at matData.ndof2 + 1 + + sum_nDOFs2 = sum_nDOFs2 + mod_nDOFs; + else: + if indx_start < mod_nDOFs: + StateOrderingIndx[indx_start:(indx_start+mod_nDOFs)] = sum_nDOFs1 + matData['NumStates2'] + np.arange(0,mod_nDOFs); # q1 starts at matData.NumStates2 + 1 + + sum_nDOFs1 = sum_nDOFs1 + mod_nDOFs; + + # reset for a new module (or new 1st-order states in the same module) + mod_nDOFs = 0; + + indx_start = i; #start of this module + lastModName = modName; + lastModOrd = matData['StateDerivOrder'][i]; + + mod_nDOFs = mod_nDOFs+1; + + + # repeat for the last module found: + if lastModOrd == 2: + mod_nDOFs = int(mod_nDOFs/2); + #print(mod_nDOFs,indx_start) + StateOrderingIndx[indx_start:(indx_start+mod_nDOFs)] = sum_nDOFs2 + np.arange(0,mod_nDOFs); # q2 starts at 1 + StateOrderingIndx[(indx_start+mod_nDOFs):matData['NumStates']] = sum_nDOFs2 + matData['ndof2'] + np.arange(0,mod_nDOFs); # q2_dot starts at matData.ndof2 + 1 + else: + StateOrderingIndx[indx_start:(indx_start+mod_nDOFs)] = sum_nDOFs1 + matData['NumStates2'] + np.arange(0,mod_nDOFs); # q1 starts at matData.NumStates2 + 1 + + #print(StateOrderingIndx) + return StateOrderingIndx + +def readOP(fid, n, defaultDerivOrder=2): + OP=[] + Var = {'RotatingFrame': [], 'DerivativeOrder': [], 'Description': []} + colNames=fid.readline().strip() + dummy= fid.readline().strip() + bHasDeriv= colNames.find('Derivative Order')>=0 + for i, line in enumerate(fid): + sp=line.strip().split() + if sp[1].find(',')>=0: + # Most likely this OP has three values (e.g. orientation angles) + # For now we discard the two other values + OP.append(float(sp[1][:-1])) + iRot=4 + else: + OP.append(float(sp[1])) + iRot=2 + Var['RotatingFrame'].append(sp[iRot]) + if bHasDeriv: + Var['DerivativeOrder'].append(int(sp[iRot+1])) + Var['Description'].append(' '.join(sp[iRot+2:]).strip()) + else: + Var['DerivativeOrder'].append(defaultDerivOrder) + Var['Description'].append(' '.join(sp[iRot+1:]).strip()) + if i>=n-1: + break + + tmp = dict() + tmp['x_op'] = OP + tmp['x_rotFrame'] = Var['RotatingFrame'] + tmp['x_DerivOrder'] = Var['DerivativeOrder'] + tmp['x_desc'] = Var['Description'] + + return tmp + + + +def readFASTMatrix(f): + name="" + m=0 + tmp=[] + for line in islice(f,1): + # copy matrix name to tmp_name + tmp_name=line.strip().split(':')[0] + # get matrix dimensions if matrix exists + if tmp_name != "": + m=int(line.strip().split(':')[1].split('x')[0]) + n=int(line.strip().split(':')[1].split('x')[1]) + + # copy matrix into tmp list + if m!=0: + name=tmp_name + for line in islice(f,m): + tmp.append([float(num) for num in line.split()]) + tmp=np.array(tmp) + return name,tmp + +def ReadFASTLinear(filename): + + def extractVal(lines, key): + for l in lines: + if l.find(key)>=0: + return l.split(key)[1].split()[0] + return None + + def readToMarker(fid, marker, nMax): + lines=[] + for i, line in enumerate(fid): + if i>nMax: + raise BrokenFormatError('`{}` not found in file'.format(marker)) + if line.find(marker)>=0: + break + lines.append(line.strip()) + return lines, line + + with open(filename) as f: + info = {} + data = {} + SetOfMatrices = 1 + info['name'] = os.path.splitext(os.path.basename(filename))[0] + + # --- + header, lastLine=readToMarker(f, 'Jacobians included', 30) + header.append(lastLine) + data['t'] = float(extractVal(header,'Simulation time:' )) + data['n_x'] = int(extractVal(header,'Number of continuous states:')) + data['n_xd'] = int(extractVal(header,'Number of discrete states:' )) + data['n_z'] = int(extractVal(header,'Number of constraint states:')) + data['n_u'] = int(extractVal(header,'Number of inputs:' )) + data['n_y'] = int(extractVal(header,'Number of outputs:' )) + bJac = extractVal(header,'Jacobians included in this file?') + if bJac: + SetOfMatrices = 2 + try: + data['Azimuth'] = float(extractVal(header,'Azimuth:')) + except: + data['Azimuth'] = None + try: + data['RotSpeed'] = float(extractVal(header,'Rotor Speed:')) # rad/s + except: + data['RotSpeed'] = np.nan + try: + data['WindSpeed'] = float(extractVal(header,'Wind Speed:')) + except: + data['WindSpeed'] = np.nan + + # --- Old method for reading + # header = [f.readline() for _ in range(17)] + # info['fast_version'] = header[1].strip() + # info['modules'] = header[2].strip() + # info['description'] = header[4].strip() + # ln=7; + # data = np.array([line.split() for line in f.readlines()]).astype(np.float) + # data['t']=np.float(header[ln].split(':')[1].strip().split(' ')[0]) + # data['RotSpeed']=np.float(header[ln+1].split(':')[1].strip().split(' ')[0]) + # data['Azimuth']=np.float(header[ln+2].split(':')[1].strip().split(' ')[0]) + # data['WindSpeed']=np.float(header[ln+3].split(':')[1].strip().split(' ')[0]) + # data['n_x']=np.float(header[ln+4].split(':')[1].strip().split(' ')[0]) + # data['n_xd']=np.float(header[ln+5].split(':')[1].strip().split(' ')[0]) + # data['n_z']=np.float(header[ln+6].split(':')[1].strip().split(' ')[0]) + # data['n_u']=np.float(header[ln+7].split(':')[1].strip().split(' ')[0]) + # data['n_y']=np.float(header[ln+8].split(':')[1].strip().split(' ')[0]) + # if header[ln+9].split('?')[1].strip()=='Yes': + # SetOfMatrices=2 + + data['Azimuth']=np.mod(data['Azimuth'],2.0*np.pi) + try: + # skip next three lines + for line in islice(f,2): + pass + if data['n_x'] > 0: + temp = readOP(f, data['n_x']) + data['x_op']=temp['x_op'] + data['x_rotFrame']=temp['x_rotFrame'] + data['x_DerivOrder']=temp['x_DerivOrder'] + data['x_desc']=temp['x_desc'] + + # skip next three lines + for line in islice(f,2): + pass + + temp = readOP(f, data['n_x'], defaultDerivOrder=2) + data['xdot_op']=temp['x_op'] + data['xdot_desc']=temp['x_desc'] + + #(number of second-order states) + data['n_x2'] = sum(1 for i in data['x_DerivOrder'] if i == 2) + else: + data['n_x2'] = 0; + + + if data['n_xd'] > 0: + # skip next three lines + for line in islice(f,2): + pass + temp = readOP(f, data['n_xd'], defaultDerivOrder=2) + data['xd_op']=temp['x_op'] + data['xd_desc']=temp['x_desc'] + if data['n_z'] > 0: + # skip next three lines + for line in islice(f,2): + pass + temp = readOP(f, data['n_z'], defaultDerivOrder=0) + data['z_op']=temp['x_op'] + data['z_desc']=temp['x_desc'] + if data['n_u'] > 0: + # skip next three lines + for line in islice(f,2): + pass + temp = readOP(f, data['n_u'], defaultDerivOrder=0) + data['u_op']=temp['x_op'] + data['u_desc']=temp['x_desc'] + data['u_rotFrame']=temp['x_rotFrame'] + if data['n_y'] > 0: + # skip next three lines + for line in islice(f,2): + pass + temp = readOP(f, data['n_y'], defaultDerivOrder=0) + data['y_op']=temp['x_op'] + data['y_desc']=temp['x_desc'] + data['y_rotFrame']=temp['x_rotFrame'] + + # skip next one line + for line in islice(f,4): + pass + + mat_names=[] + while True: + name,mat=readFASTMatrix(f) + if not name: + break; + mat_names.append(name) + data[name]=mat + + return data, info + + except (ValueError, AssertionError): + raise + + +def get_Mats(FileNames, verbose=True, removeTwrAzimuth=False): + NAzimStep = len(FileNames) + data = [None]*NAzimStep + # --- Read all files + for iFile, filename in enumerate(FileNames): + data[iFile],_= ReadFASTLinear(FileNames[iFile]); + Azimuth = np.array([d['Azimuth'] for d in data])*180/np.pi + # --- Sort by azimuth, not required, but filenames are not sorted, nor are the lin times azimuth + ISort = np.argsort(Azimuth) + Azimuth = Azimuth[ISort] + data = [data[i] for i in ISort] + FileNames = [FileNames[i] for i in ISort] + if removeTwrAzimuth: + IDiscard = [i for i,a in enumerate(Azimuth) if np.any(np.abs(np.array([60,180,300])-a)<4) ] + n=len(FileNames[0]); sfmt='{:'+str(n+2)+'s}' + for i in IDiscard: + print(' discard: '+sfmt.format(FileNames[i]) + ' (psi={:5.1f})'.format(Azimuth[i])) + data = [d for i,d in enumerate(data) if i not in IDiscard] + NAzimStep= len(data) + + + matData={} + matData['NAzimStep']= NAzimStep; + + # --- Using last azimuth to initial some of the "matData" variables + # Input data from linearization files + matData['NumStates'] = int(data[NAzimStep-1]['n_x']); + matData['NumStates2'] = int(data[NAzimStep-1]['n_x2']); + + #return data, matData + matData['ndof1'] = int(matData['NumStates'] - matData['NumStates2']); # number of first-order states = number of first-order DOFs + matData['ndof2'] = int(data[NAzimStep-1]['n_x2'] / 2); # half the number of second-order states = number of second-order DOFs + #% matData['ndof'] = matData['ndof2'] + matData['ndof1']; #half the number of second-order states plus the number of first-order states (i.e., states that aren't derivatives) + + matData['NumInputs'] = int(data[NAzimStep-1]['n_u']); + matData['NumOutputs'] = int(data[NAzimStep-1]['n_y']); + + # allocate space for these variables + matData['Azimuth'] = np.zeros(NAzimStep); + matData['Omega'] = np.zeros(NAzimStep); + matData['OmegaDot'] = np.zeros(NAzimStep); + + matData['WindSpeed'] = np.zeros(NAzimStep)*np.nan; + + if matData['NumStates'] > 0: + matData['DescStates'] = data[NAzimStep-1]['x_desc']; + matData['StateDerivOrder'] = data[NAzimStep-1]['x_DerivOrder']; + matData['xdop'] = np.zeros((matData['NumStates'], NAzimStep)) + matData['xop'] = np.zeros((matData['NumStates'], NAzimStep)) + matData['A'] = np.zeros((matData['NumStates'], matData['NumStates'], NAzimStep)) + + if matData['NumInputs'] > 0: + matData['DescCntrlInpt'] = data[NAzimStep-1]['u_desc']; + matData['u_op'] = np.zeros((matData['NumInputs'],NAzimStep)) + if matData['NumStates']>0: + matData['B'] = np.zeros((matData['NumStates'], matData['NumInputs'],NAzimStep)) + + if matData['NumOutputs'] > 0: + matData['DescOutput'] = data[NAzimStep-1]['y_desc'] + matData['y_op'] = np.zeros((matData['NumOutputs'],NAzimStep)) + + if matData['NumStates'] > 0: + matData['C'] = np.zeros((matData['NumOutputs'], matData['NumStates'], NAzimStep)) + if matData['NumInputs'] > 0: + matData['D'] = np.zeros((matData['NumOutputs'], matData['NumInputs'], NAzimStep)) + + # Reorder state matrices so that they follow the {q2, q2_dot, q1} + # format that is assumed in the MBC3 equations. + if matData['NumStates'] > 0: + # keep StateOrderingIndx for applying inverse of MBC3 later + # (to visualize mode shapes) + matData['StateOrderingIndx'] = getStateOrderingIndx(matData) + + sortedIndx=matData['StateOrderingIndx'] + #print(sortedIndx) + x_rotFrame= reOrderByIdx_1D(data[NAzimStep-1]['x_rotFrame'],sortedIndx) + matData['DescStates'] = reOrderByIdx_1D(data[NAzimStep-1]['x_desc'],sortedIndx) + matData['StateDerivOrder'] = reOrderByIdx_1D(data[NAzimStep-1]['x_DerivOrder'],sortedIndx) + + # --- Store file data into matData + for iFile in np.arange(0,NAzimStep): + matData['Omega'][iFile] = data[iFile]['RotSpeed'] ; + matData['Azimuth'][iFile] = data[iFile]['Azimuth']*180/np.pi; + + matData['WindSpeed'][iFile] = data[iFile]['WindSpeed'] + + if 'A' in data[iFile]: + matData['A'][:,:,iFile]=reOrderByIdx_2D(data[iFile]['A'],sortedIndx,sortedIndx) + #print('size of matData[A] for file ', iFile, ' is :', matData['A'].shape) + if 'B' in data[iFile]: + matData['B'][:,:,iFile]=reOrderByIdx_2D(data[iFile]['B'],sortedIndx,np.arange(0,matData['NumStates'])) + #print('size of matData[B] for file ', iFile, ' is :', matData['B'].shape) + if 'C' in data[iFile]: + matData['C'][:,:,iFile]=reOrderByIdx_2D(data[iFile]['C'],np.arange(0,matData['NumStates']),sortedIndx) + #print('size of matData[C] for file ', iFile, ' is :', matData['C'].shape) + if 'D' in data[iFile]: + matData['D'][:,:,iFile]=reOrderByIdx_2D(data[iFile]['D'],sortedIndx,sortedIndx) + #print('size of matData[D] for file ', iFile, ' is :', matData['D'].shape) + + if 'x_op' in data[iFile]: + matData['xop'][:,iFile] = reOrderByIdx_1D(data[iFile]['x_op'],sortedIndx) + #print('matData[xop] for file ', iFile,' is :',(matData['xop'][:,iFile])) + if 'xdot_op' in data[iFile]: + matData['xdop'][:,iFile] = reOrderByIdx_1D(data[iFile]['xdot_op'],sortedIndx) + #print('matData[xdop] for file ', iFile,' is :',(matData['xdop'][:,iFile])) + + if 'u_op' in data[iFile]: + matData['u_op'][:,iFile] = data[iFile]['u_op'] + + if 'y_op' in data[iFile]: + matData['y_op'][:,iFile] = data[iFile]['y_op'] + + + # Find the azimuth-averaged linearized 1st order state matrices: + if 'A' in matData: + matData['Avgxdop'] = np.mean(matData['xdop'],axis=1) + matData['Avgxop'] = np.mean(matData['xop'], axis=1) + matData['AvgA'] = np.mean(matData['A'],axis=2) + + #print(matData['AvgA']) + if 'B' in matData: + matData['AvgB'] = np.mean(matData['B'],axis=2) + + if 'C' in matData: + matData['AvgC'] = np.mean(matData['C'],axis=2) + + if 'D' in matData: + matData['AvgD'] = np.mean(matData['D'],axis=2) + + + foundED = True; + for i in range(matData['ndof2']): + # find the starting index of the string 'DOF_GeAz' + if (matData['DescStates'][i].find('DOF_GeAz') != -1): + matData['Omega'] = matData['xdop'][i,:] + matData['OmegaDot'] = matData['xdop'][i+matData['ndof2'],:] + foundED = True + break + + for i in range(matData['ndof2']): + # find the starting index of the string 'DOF_DrTr' + if (matData['DescStates'][i].find('DOF_DrTr') != -1): + matData['Omega'] = matData['Omega'] + matData['xdop'][i,:] #This always comes after DOF_GeAz so let's just add it here (it won't get written over later). + matData['OmegaDot'] = matData['OmegaDot'] + matData['xdop'][i+matData['ndof2'],:] + foundED = True + break + + if not foundED: + for i in range(matData['ndof2']): + # find the starting index of the string 'Gearbox_Rot' + if (matData['DescStates'][i].find('MBD Gearbox_Rot') != -1): + matData['Omega'] = matData['xdop'][i,:] + matData['OmegaDot'] = matData['xdop'][i+matData['ndof2'],:] + break + + #print("\n".join(matData['DescStates'])) + #exit() + # ----------- Find multi-blade coordinate (MBC) transformation indices ---- + + # Find the indices for, state triplets in the rotating frame + # (note that we avoid the "first time derivative" states) + if matData['ndof2'] > 0: + matData['RotTripletIndicesStates2'], matData['n_RotTripletStates2'] = findBladeTriplets(x_rotFrame[0:matData['ndof2']],matData['DescStates'][0:matData['ndof2']], verbose=verbose) + else: + matData['RotTripletIndicesStates2'] = []; + matData['n_RotTripletStates2'] = 0; + + if matData['ndof1'] > 0: + matData['RotTripletIndicesStates1'], matData['n_RotTripletStates1'] = findBladeTriplets( x_rotFrame[matData['NumStates2']:] ,matData['DescStates'][matData['NumStates2']:] , verbose=verbose); + else: + matData['RotTripletIndicesStates1'] = []; + matData['n_RotTripletStates1'] = 0; + + # Find the indices for control input triplets in the rotating frame: + if matData['NumInputs'] > 0: + matData['RotTripletIndicesCntrlInpt'], matData['n_RotTripletInputs'] = findBladeTriplets(data[0]['u_rotFrame'],matData['DescCntrlInpt'], verbose=verbose ); + else: + matData['RotTripletIndicesCntrlInpt'] = []; + matData['n_RotTripletInputs'] = 0; + + # Find the indices for output measurement triplets in the rotating frame: + if (matData['NumOutputs'] > 0 ): + matData['RotTripletIndicesOutput'], matData['n_RotTripletOutputs'] = findBladeTriplets(data[0]['y_rotFrame'],matData['DescOutput'], verbose=verbose ); + else: + matData['RotTripletIndicesOutput'] = []; + matData['n_RotTripletOutputs'] = 0; + + return matData, data + +#matData,FAST_linData=get_Mats(FileNames) + +#print("\n".join(matData['DescStates'])) + +# print(info) +# print(data['t']) +# print(data['RotSpeed']) +# print(data['Azimuth']) +# print(data['WindSpeed']) +# print(data['n_x']) +# print(data['n_xd']) +# print(data['n_z']) +# print(data['n_u']) +# print(data['n_y']) +# print(data['n_x2']) +# print(data['x_op']) +# print(data['x_desc']) +# print(data['x_rotFrame']) +# print(data['xdot_op']) +# print(data['xdot_desc']) +# print(data['u_op']) +# print(data['u_desc']) +# print(data['u_rotFrame']) +# print(data['y_op']) +# print(data['y_desc']) +# print(data['y_rotFrame']) diff --git a/rosco/toolbox/linear/linear_models.py b/rosco/toolbox/linear/linear_models.py index b1ebfb78..f2004a7a 100644 --- a/rosco/toolbox/linear/linear_models.py +++ b/rosco/toolbox/linear/linear_models.py @@ -11,15 +11,9 @@ import control as co import matplotlib.pyplot as plt import multiprocessing as mp -from itertools import chain from scipy.io import loadmat -try: - import pyFAST.linearization.mbc as mbc -except ImportError: - import weis.control.mbc.mbc3 as mbc -except ImportError: - raise ImportError('Unable to load mbc3 from pyFAST or WEIS') +import rosco.toolbox.linear.mbc3 as mbc class LinearTurbineModel(object): @@ -710,7 +704,7 @@ def run_mbc3(fnames): Helper function to run mbc3 ''' print('Loading linearizations from:', ''.join(fnames[0].split('.')[:-2])) - MBC, matData = mbc.fx_mbc3(fnames, verbose=False) + MBC, matData, _ = mbc.fx_mbc3(fnames, verbose=False) return MBC, matData diff --git a/rosco/toolbox/linear/mbc3.py b/rosco/toolbox/linear/mbc3.py new file mode 100644 index 00000000..e80956e8 --- /dev/null +++ b/rosco/toolbox/linear/mbc3.py @@ -0,0 +1,1071 @@ +############################################################# +# Perform MBC3 Analysis and plot Campbell Diagram # +############################################################# + +import rosco.toolbox.linear.getMats as gm +import numpy as np +import scipy.linalg as scp +import re +import os +import sys +# import plotCampbellData as pCD + +def getScaleFactors(DescStates, TowerLen, BladeLen): + + ScalingFactor = np.ones(len(DescStates)) + + # look at the state description strings for tower and blade + # translational dofs: + for i in range(len(ScalingFactor)): + + # look for blade dofs: + if DescStates[i].find('blade')!=-1 or DescStates[i].find('Blade')!=-1: + if DescStates[i].find('rotational')==-1: # make sure this isn't a rotational dof from BeamDyn + ScalingFactor[i] = 1.0/BladeLen; + #print(DescStates[i]) + + # look for tower translational dofs: + if DescStates[i].find('tower')!=-1 or DescStates[i].find('Tower')!=-1: + ScalingFactor[i] = 1.0/TowerLen; + + # look for blade dofs: + elif DescStates[i].find('blade')!=-1 or DescStates[i].find('Blade')!=-1: + if DescStates[i].find('rotational')==-1: # make sure this isn't a rotational dof from BeamDyn + ScalingFactor[i] = 1.0/BladeLen; + + return ScalingFactor + + +def IdentifyModes(CampbellData): + """ + Attempts to perform an identification of the modes. + For now, the method is based on the energy content of the modes, and the state descriptions where the energy is maximum + + Original contribution by: Srinivasa B. Ramisett, ramisettisrinivas@yahoo.com, http://ramisetti.github.io + """ + + # --- Looking at states descriptions (of first run, first mode), to see if we are offshore. + # NOTE: DescStates is likely the same for all modes + DescStates = CampbellData[0]['Modes'][0]['DescStates'] + hasHeave = any(['heave' in s.lower() for s in DescStates]) + hasSurge = any(['surge' in s.lower() for s in DescStates]) + hasSway = any(['sway' in s.lower() for s in DescStates]) + hasYaw = any(['platform yaw' in s.lower() for s in DescStates]) + hasRoll = any(['platform roll tilt' in s.lower() for s in DescStates]) + hasPitch = any(['platform pitch tilt' in s.lower() for s in DescStates]) + hasEdge1Col = any(['1st edgewise bending-mode dof of blade collective' in s.lower() for s in DescStates]) + + # --- Setting up a list of modes with + modesDesc = [] + modesDesc.append( ['Generator DOF (not shown)' , 'ED Variable speed generator DOF, rad'] ) + # Platform DOFs + if hasSurge: + modesDesc.append( ['Platform surge', 'ED Platform horizontal surge translation DOF, m'] ) + if hasSway: + modesDesc.append( ['Platform sway', 'ED Platform horizontal sway translation DOF, m'] ) + if hasHeave: + modesDesc.append( ['Platform heave', 'ED Platform vertical heave translation DOF, m'] ) + if hasRoll: + modesDesc.append( ['Platform roll', 'ED Platform roll tilt rotation DOF, rad'] ) + if hasPitch: + modesDesc.append( ['Platform pitch', 'ED Platform pitch tilt rotation DOF, rad'] ) + if hasYaw: + modesDesc.append( ['Platform yaw', 'ED Platform yaw rotation DOF, rad'] ) + + modesDesc.append( ['1st Tower FA' , 'ED 1st tower fore-aft bending mode DOF, m'] ) + modesDesc.append( ['1st Tower SS' , 'ED 1st tower side-to-side bending mode DOF, m'] ) + modesDesc.append( ['1st Blade Flap (Regressive)' , 'ED 1st flapwise bending-mode DOF of blade (sine|cosine), m', r'Blade (sine|cosine) finite element node \d rotational displacement in Y, rad'] ) + modesDesc.append( ['1st Blade Flap (Collective)' , 'ED 1st flapwise bending-mode DOF of blade collective, m', r'Blade collective finite element node \d rotational displacement in Y, rad'] ) + modesDesc.append( ['1st Blade Flap (Progressive)' , 'ED 1st flapwise bending-mode DOF of blade (sine|cosine), m'] ) # , ... # 'Blade (sine|cosine) finite element node \d rotational displacement in Y, rad'] + modesDesc.append( ['1st Blade Edge (Regressive)' , 'ED 1st edgewise bending-mode DOF of blade (sine|cosine), m', r'Blade (sine|cosine) finite element node \d rotational displacement in X, rad'] ) + if hasEdge1Col: + modesDesc.append(['1st Blade Edge (Collective)', 'ED 1st edgewise bending-mode DOF of blade collective, m'] ) # , ... # 'Blade (sine|cosine) finite element node \d rotational displacement in Y, rad'] + modesDesc.append( ['1st Blade Edge (Progressive)' , 'ED 1st edgewise bending-mode DOF of blade (sine|cosine), m'] ) + modesDesc.append( ['1st Drivetrain Torsion' , 'ED Drivetrain rotational-flexibility DOF, rad'] ) + modesDesc.append( ['2nd Tower FA' , 'ED 2nd tower fore-aft bending mode DOF, m'] ) + modesDesc.append( ['2nd Tower SS' , 'ED 2nd tower side-to-side bending mode DOF, m'] ) + modesDesc.append( ['2nd Blade Flap (Regressive)' , 'ED 2nd flapwise bending-mode DOF of blade (sine|cosine), m'] ) + modesDesc.append( ['2nd Blade Flap (Collective)' , 'ED 2nd flapwise bending-mode DOF of blade collective, m', r'Blade collective finite element node \d rotational displacement in Y, rad'] ) + modesDesc.append( ['2nd Blade Flap (Progressive)' , 'ED 2nd flapwise bending-mode DOF of blade (sine|cosine), m'] ) + modesDesc.append( ['Nacelle Yaw (not shown)' , 'ED Nacelle yaw DOF, rad'] ) + + + nModes = int(len(modesDesc)) + nRuns = int(len(CampbellData)) + modeID_table=np.zeros((nModes,nRuns)).astype(int) + + + + def doesDescriptionMatch(description, listOfModePatterns): + """ loop through all mode desrption """ + for iModePattern, modeIDdescList in enumerate(listOfModePatterns): + modeIDName = modeIDdescList[0] + patternList = modeIDdescList[1:] # list of patterns for a given mode + found, pattern = doesDescriptionMatchPatternList(description, patternList) + if found: + return True, iModePattern, pattern + return False, -1, '' + + def doesDescriptionMatchPatternList(description, patternList): + """ loop through all patterns to find a match """ + for pattern in patternList: + # Looking for targetDesc into description + if re.search(pattern ,description, re.IGNORECASE)!=None: + return True, pattern + return False, '' + + + verbose=False + Levels=[1,2,3] + + + # --- Loop on operating points/Runs + for i in range(nRuns): + Modes = CampbellData[i]['Modes'] + nModes = len(Modes) + # Array of logical, False for Modes that are not identified + modesIdentified = [False] * nModes + modesSkipped = [False] * nModes + #verbose=verbose and i==0 # only display for first mode for now.. + #verbose=i==1 + + # --- Give an index to each mode so that we can easily identify them + for im, mode in enumerate(Modes): + mode['index']=im + + # --- Skip modes based on simple criteria + for im, mode in enumerate(Modes): + if mode['NaturalFreq_Hz'] < 1e-5 or mode['DampingRatio'] > 0.98: + modesSkipped[im]=True + + + + modesDescLocal = modesDesc.copy() + + # --- Level 1 - Find well-defined modes (modes which have only one Max) + if 1 in Levels: + for im, mode in enumerate(Modes): + if modesIdentified[im] or modesSkipped[im]: + continue # skip this mode since it has already been identified + stateMax=np.argwhere((mode['StateHasMaxAtThisMode']==1)).flatten() + if len(stateMax)==1: + description = mode['DescStates'][stateMax[0]] + if description.startswith('AD'): + # We skipp the pure "AD" modes + modesSkipped[im] = True + continue + found, modeID, patternMatched = doesDescriptionMatch(description, modesDescLocal) + if found and modeID_table[modeID,i]==0: + modesDescLocal[modeID] = [None,] # we empty this mode patternlist so that it cannot be matched again + modesIdentified[im] = True + modeID_table[modeID,i] = im+1 # Using Matlab Indexing + if verbose: + print('L1 Mode {} identified using pattern: {}'.format(im+1,patternMatched)) + print(' String was: ', description) + #else: + # print('>>> Cannot identify mode with description {}. Update MBC3 script.'.format(description)) + + # --- Level 2 - Find modes with several max - Looping on mode pattern to respect expected frequency order + if 2 in Levels: + for modeID, modeIDdescList in enumerate(modesDescLocal): + modeIDName = modeIDdescList[0] + patternList = modeIDdescList[1:] + # Skip modes already identified above + if modeID_table[modeID,i]>0: + continue + if verbose: + print('------------------------- LEVEL 2 - LOOKING FOR MODE ',modeIDName) + + found = False; + # --- Loop on all non-identified modes in increasing order + im = 0; + for im, mode in enumerate(Modes): + if modesIdentified[im] or modesSkipped[im]: + continue # move to next mode + # List of component descriptions where this mode has maximum values + stateMax=np.argwhere((mode['StateHasMaxAtThisMode']==1)).flatten() + descriptions = np.array(mode['DescStates'])[stateMax] + descriptions = descriptions[:7] # we keep only the first 7 descriptions + descriptionsED = [d for d in descriptions if d.startswith('ED')] + descriptionsBD = [d for d in descriptions if d.startswith('BD')] + descriptionsAD = [d for d in descriptions if d.startswith('AD')] + descriptionsMisc = [d for d in descriptions if d not in descriptionsED+descriptionsBD+descriptionsAD] + descriptions = descriptionsED+descriptionsBD+descriptionsMisc # NOTE: we skipp AD descriptions + j = 0; + for description in descriptions: + found, pattern = doesDescriptionMatchPatternList(description, patternList) + if found: + if verbose: + print('L2 Mode {} identified using pattern {}'.format(im+1,pattern)) + modeID_table[modeID,i] = im+1 # Using Matlab Indexing + modesDescLocal[modeID] = [None,] # we empty this mode patternlist so that it cannot be matched again + modesIdentified[im] = True; + break + if found: + break + if verbose: + print('>> modeIDTable',modeID_table[:,i]) + # We disqualify modes that had max and that didn't match anything: + for im, mode in enumerate(Modes): + if modesIdentified[im] or modesSkipped[im]: + continue # move to next mode + stateMax=np.argwhere((mode['StateHasMaxAtThisMode']==1)).flatten() + if len(stateMax)>=1: + modesSkipped[im]=True + shortdescr = CampbellData[i]['ShortModeDescr'][im] + if shortdescr.find('ED')>=0: + print('>>>> short', CampbellData[i]['ShortModeDescr'][im]) + print('>>>> Problem in IndentifyModes. ED DOF found in level 2') + # import pdb; pdb.set_trace() + + if 3 in Levels: + # --- Level 3 - Try our best for modes with no max + # Loop on modes to be identified + for modeID, modeIDdescList in enumerate(modesDescLocal): + modeIDName = modeIDdescList[0] + patternList = modeIDdescList[1:] + + # Skip modes already identified above + if modeID_table[modeID,i]>0: + continue + if verbose: + print('------------------------- LEVEL 3 - LOOKING FOR MODE ',modeIDName) + + found = False; + # --- Loop on all non-identified modes in increasing order + im = 0; + while not found and im < nModes: # Loop on modes + mode = Modes[im] + if modesIdentified[im] or modesSkipped[im]: + pass + else: + # --- Otherwise, use as mode descriptions the other ones. Seems weird + stateMax=np.argwhere((mode['StateHasMaxAtThisMode']==0)).flatten() + descriptions = np.array(mode['DescStates'])[stateMax] + ADcounts = np.sum([s.startswith('AD') for s in descriptions[:5]]) + + descriptions2= np.array(mode['DescStates'])[mode['StateHasMaxAtThisMode']] + if len(descriptions2) == 0: + noMax=True + descriptions3 = mode['DescStates'][:5] + else: + noMax=False +# import pdb; pdb.set_trace() + if ADcounts<5: + descriptions=[d for d in descriptions if not d.startswith('AD')] +# descriptionsED = [d for d in descriptions if d.startswith('ED')] +# descriptionsBD = [d for d in descriptions if d.startswith('BD')] +# descriptionsAD = [d for d in descriptions if d.startswith('AD')] +# descriptionsMisc = [d for d in descriptions if d not in descriptionsED+descriptionsBD+descriptionsAD] +# descriptions = descriptionsED+descriptionsBD+descriptionsMisc # NOTE: we skipp AD descriptions + descriptions = descriptions[:5] # we keep only the first 7 descriptions + if verbose: + print('>>> Mode',mode['index'], modesIdentified[im], modesSkipped[im]) + print('>>>> descr', [replaceModeDescription(s) for s in descriptions]) + print('>>>> short', CampbellData[i]['ShortModeDescr'][im]) + else: + descriptions=[] +# #descriptions = descriptions[:7] # we keep only the first 7 descriptions + + j = 0; + while not found and j < len(descriptions): + j = j + 1; + if not found: + for targetDesc in patternList: + # Looking for targetDesc into list of descriptions + if re.search(targetDesc ,descriptions[j-1],re.IGNORECASE)!=None: + modeID_table[modeID,i] = im+1 # Using Matlab Indexing + if verbose: + print('L3 Mode {} identified as {}'.format(im+1,targetDesc)) + print(' String was: ', descriptions[j-1]) + modesIdentified[im] = True; + found = True; + break; + im=im+1; # increment counter + if verbose: + print('>> modeIDTable',modeID_table[:,i]) + + if verbose: + print('---------- Summary') + for j in np.arange(len(modeID_table)): + print('{:32s} {:d}'.format(modesDesc[j][0],modeID_table[j,i])) + print('---------- ') + + + return modeID_table,modesDesc + +def IdentifiedModesDict(CampbellData,modeID_table,modesDesc): + """ + To be called with the results of IdentifyModes. + Create a list of dictionaries to more easily interprete the result + """ + nOP = modeID_table.shape[1] + modesInfoPerOP=[] + for iOP in np.arange(nOP): + modesInfo={} + for i in np.arange(len(modesDesc)): + desc = modesDesc[i][0] + ID = int(modeID_table[i,iOP]) + if ID==0: + freq=np.nan + damp=np.nan + cont='' + else: + freq = np.around(CampbellData[iOP]['Modes'][ID-1]['NaturalFreq_Hz'],5) + damp = np.around(CampbellData[iOP]['Modes'][ID-1]['DampingRatio'],5) + cont = CampbellData[iOP]['ShortModeDescr'][ID-1] + modesInfo[desc]={'ID':ID,'f0':freq,'zeta':damp,'cont':cont} + modesInfoPerOP.append(modesInfo) + return modesInfoPerOP + +def campbell_diagram_data(mbc_data, BladeLen, TowerLen): + """ + + Original contribution by: Srinivasa B. Ramisett, ramisettisrinivas@yahoo.com, http://ramisetti.github.io + """ + + CampbellData={} + usePercent = False; + # + # mbc_data.eigSol = eiganalysis(mbc_data.AvgA); + ndof = mbc_data['ndof2'] + mbc_data['ndof1']; #size(mbc_data.AvgA,1)/2; # number of translational states + nModes = len(mbc_data['eigSol']['Evals']) + #print(nModes) + DescStates = PrettyStateDescriptions(mbc_data['DescStates'], mbc_data['ndof2'], mbc_data['performedTransformation']); + + ## store indices of max mode for state and to order natural frequencies + #StatesMaxMode_vals = np.amax(mbc_data['eigSol']['MagnitudeModes'],axis=1); # find which mode has the maximum value for each state (max of each row before scaling) + StatesMaxMode = np.argmax(mbc_data['eigSol']['MagnitudeModes'],axis=1); # find which mode has the maximum value for each state (max of each row before scaling) + SortedFreqIndx = np.argsort((mbc_data['eigSol']['NaturalFreqs_Hz']).flatten(),kind="heapsort"); + #print(SortedFreqIndx) + + + if BladeLen!=0 or TowerLen!=0: + ## get the scaling factors for the mode rows + ScalingFactor = getScaleFactors(DescStates, TowerLen, BladeLen); + + ## scale the magnitude of the modes by ScalingFactor (for consistent units) + # and then scale the columns so that their maximum is 1 + + ModesMagnitude = np.matmul(np.diag(ScalingFactor), mbc_data['eigSol']['MagnitudeModes']); # scale the rows + #print(ModesMagnitude) + + CampbellData['ScalingFactor'] = ScalingFactor; + else: + ModesMagnitude = mbc_data['eigSol']['MagnitudeModes']; + + if usePercent: + scaleCol = np.sum( ModesMagnitude )/100; # find the sum of the column, and multiply by 100 (divide here) to get a percentage + else: + scaleCol = np.amax(ModesMagnitude,axis=0); #find the maximum value in the column, so the first element has value of 1 + + ModesMagnitude = np.matmul(ModesMagnitude,np.diag(1./scaleCol)) # scale the columns + + # --- Summary data (array for all modes) + CampbellData['NaturalFreq_Hz'] = mbc_data['eigSol']['NaturalFreqs_Hz'][SortedFreqIndx] + CampbellData['DampingRatio'] = mbc_data['eigSol']['DampRatios'][SortedFreqIndx] + CampbellData['RotSpeed_rpm'] = mbc_data['RotSpeed_rpm'] + if 'WindSpeed' in mbc_data: + CampbellData['WindSpeed'] = mbc_data['WindSpeed'] + + #print(ModesMagnitude) + + # --- Storing data per mode + CampbellData['Modes']=[] + + for i in range(nModes): + CData={} + CData['NaturalFreq_Hz'] = mbc_data['eigSol']['NaturalFreqs_Hz'][SortedFreqIndx[i]][0] + CData['DampedFreq_Hz'] = mbc_data['eigSol']['DampedFreqs_Hz'][SortedFreqIndx[i]][0]; + CData['DampingRatio'] = mbc_data['eigSol']['DampRatios'][SortedFreqIndx[i]][0]; + + + #print(np.argsort(ModesMagnitude[:,SortedFreqIndx[0]])[::-1]) + # sort indices in descending order + sort_state = np.argsort( ModesMagnitude[:,SortedFreqIndx[i]])[::-1]; + + #print(type(sort_state)) + CData['DescStates']=[DescStates[i] for i in sort_state] + CData['MagnitudePhase']=ModesMagnitude[sort_state,SortedFreqIndx[i]]; + Phase = mbc_data['eigSol']['PhaseModes_deg'][sort_state,SortedFreqIndx[i]]; + # if the phase is more than +/- 90 degrees different than the first + # one (whose value == 1 or is the largest %), we'll stick a negative value on the magnitude: + Phase = np.mod(Phase, 360); + + CData['PhaseDiff'] = np.mod( Phase - Phase[0], 360); # difference in range [0, 360) + PhaseIndx = CData['PhaseDiff'] > 180; + CData['PhaseDiff'][PhaseIndx] = CData['PhaseDiff'][PhaseIndx] - 360; # move to range (-180, 180] + + if ~usePercent: + PhaseIndx = CData['PhaseDiff'] > 90; + CData['MagnitudePhase'][PhaseIndx] = -1*CData['MagnitudePhase'][PhaseIndx]; + CData['PhaseDiff'][PhaseIndx] = CData['PhaseDiff'][PhaseIndx] - 180; + + PhaseIndx = CData['PhaseDiff'] <= -90; + CData['MagnitudePhase'][PhaseIndx] = -1*CData['MagnitudePhase'][PhaseIndx]; + CData['PhaseDiff'][PhaseIndx] = CData['PhaseDiff'][PhaseIndx] + 180; + + #print(CData['MagnitudePhase']) + #print(CData['PhaseDiff']) + + CData['StateHasMaxAtThisMode'] = np.ones(ndof, dtype=bool); + ix = (StatesMaxMode == SortedFreqIndx[i]); + tmp=ix[sort_state] + CData['StateHasMaxAtThisMode']=tmp + + #print(CData['StateHasMaxAtThisMode']) + #print(CData['NaturalFreq_Hz']) + CampbellData['Modes'].append(CData) + + #print(CampbellData[0]['MagnitudePhase']) + # Adding short description to summary + CampbellData['ShortModeDescr'] =[extractShortModeDescription(CampbellData['Modes'][i]) for i in range(nModes)] + + + CampbellData['nColsPerMode'] = 5; + CampbellData['ModesTable'] = {} + + for i in range(nModes): + colStart = i*CampbellData['nColsPerMode']; + CampbellData['ModesTable'][1, colStart+1 ] = 'Mode number:'; + CampbellData['ModesTable'][1, colStart+2 ] = i; + + CampbellData['ModesTable'][2, colStart+1 ] = 'Natural (undamped) frequency (Hz):'; + CampbellData['ModesTable'][2, colStart+2 ] = CampbellData['Modes'][i]['NaturalFreq_Hz'] + + CampbellData['ModesTable'][3, colStart+1 ] = 'Damped frequency (Hz):'; + CampbellData['ModesTable'][3, colStart+2 ] = CampbellData['Modes'][i]['DampedFreq_Hz'] + + CampbellData['ModesTable'][4, colStart+1 ] = 'Damping ratio (-):'; + CampbellData['ModesTable'][4, colStart+2 ] = CampbellData['Modes'][i]['DampingRatio'] + + CampbellData['ModesTable'][5, colStart+1 ] = 'Mode ' + str(i) + ' state description'; + CampbellData['ModesTable'][5, colStart+2 ] = 'State has max at mode ' + str(i); + if usePercent: + CampbellData['ModesTable'][5, colStart+3 ] = 'Mode ' + str(i) + ' contribution (%)'; + else: + CampbellData['ModesTable'][5, colStart+3 ] = 'Mode ' + str(i) + ' signed magnitude'; + + CampbellData['ModesTable'][5, colStart+4 ] = 'Mode ' + str(i) + ' phase (deg)'; + + # need to cross check these 4 lines + CampbellData['ModesTable'][6,colStart+1] = CampbellData['Modes'][i]['DescStates']; + CampbellData['ModesTable'][6,colStart+2] = CampbellData['Modes'][i]['StateHasMaxAtThisMode']; + CampbellData['ModesTable'][6,colStart+3] = CampbellData['Modes'][i]['MagnitudePhase']; + CampbellData['ModesTable'][6,colStart+4] = CampbellData['Modes'][i]['PhaseDiff']; + + #print(CampbellData['ModesTable']) + return CampbellData + + + +def extractShortModeDescription(mode): + """ + Look at which modes have max, append these description, perform shortening substitution + The description starts with noMax if no maximum exsits in this mode + The ElastoDyn modes are placed first + """ + descriptions = np.array(mode['DescStates'])[mode['StateHasMaxAtThisMode']] + if len(descriptions) == 0: + noMax=True + descriptions = mode['DescStates'][:5] + else: + noMax=False + sED = [s for s in descriptions if s.startswith('ED')] + sBD = [s for s in descriptions if s.startswith('BD')] + sMisc = [s for s in descriptions if s not in sED+sBD] + sAll = [replaceModeDescription(s) for s in sED+sBD+sMisc] + shortdescr = ' - '.join(sAll) + if noMax: + shortdescr = 'NoMax - ' + shortdescr + return shortdescr + + +def replaceModeDescription(s): + """ Perform substitutions to make the mode description shorter""" + s = s.replace('Blade','Bld') + s = s.replace('blade','Bld') + s = s.replace('First time derivative of','d/dt of') + s = s.replace('fore-aft bending mode DOF, m','FA') + s = s.replace('side-to-side bending mode DOF, m','SS') + s = s.replace('bending-mode DOF of Bld ','') + s = s.replace(' rotational-flexibility DOF, rad','-rot') + s = s.replace('rotational displacement in ','rot') + s = s.replace('translational displacement in ','trans') + s = s.replace('Platform horizontal surge translation DOF','Platform surge') + s = s.replace('Platform vertical heave translation DOF','Platform heave') + s = s.replace('Platform pitch tilt rotation DOF','Platform pitch') + s = s.replace(', rad','') + s = s.replace(', m','') + s = s.replace('finite element node ','N') + s = s.replace('cosine','cos') + s = s.replace('sine','sin') + s = s.replace('flapwise','FLAP') + s = s.replace('edgewise','EDGE') + s = s.replace('collective','coll.') + s = s.replace('rotZ','TORS-ROT') + s = s.replace('transX','FLAP-DISP') + s = s.replace('transY','EDGE-DISP') + s = s.replace('rotX','EDGE') + s = s.replace('rotY','FLAP') + return s + + +def PrettyStateDescriptions(DescStates, ndof2, performedTransformation): + idx=np.array(list(range(0,ndof2))+list(range(ndof2*2+1,len(DescStates)))) + tmpDS = [DescStates[i] for i in idx] + + if performedTransformation: + key_vals=[['BD_1','Blade collective'],['BD_2','Blade cosine'],['BD_3','Blade sine'],['blade 1','blade collective'], ['blade 2','blade cosine'], ['blade 3','Blade sine '], + ['PitchBearing1','Pitch bearing collective '], ['PitchBearing2','Pitch bearing cosine '], ['PitchBearing3','Pitch bearing sine '], + ['at Blade', 'at blade'], + ['of Blade', 'of blade'], + ] + # Replace Substrings from String List + sub = dict(key_vals) + for key, val in sub.items(): + for idx, ele in enumerate(tmpDS): + if key in ele: + tmpDS[idx] = ele.replace(key, val) + + StateDesc=tmpDS + else: + StateDesc = tmpDS + + for i in range(len( StateDesc )): + First = re.split(r'\(',StateDesc[i],2) + Last = re.split(r'\)',StateDesc[i],2) + + if len(First)>0 and len(Last)>0 and len(First[0]) != len(StateDesc[i]) and len( Last[-1] ) != len(StateDesc[i]): + StateDesc[i] = (First[0]).strip() + Last[-1]; + #print(StateDesc[i]) + + return StateDesc + + +def get_tt_inverse(sin_col, cos_col): + + c1 = cos_col[0]; + c2 = cos_col[1]; + c3 = cos_col[2]; + + s1 = sin_col[0]; + s2 = sin_col[1]; + s3 = sin_col[2]; + + + ttv = [ [c2*s3 - s2*c3, c3*s1 - s3*c1, c1*s2 - s1*c2], + [ s2 - s3 , s3 - s1, s1 - s2 ], + [ c3 - c2 , c1 - c3, c2 - c1 ] ] + ttv = ttv/(1.5*np.sqrt(3)); + + return ttv + +def get_new_seq(rot_triplet,ntot): +# rot_triplet is size n x 3 + #print('rot_triplet ', len(rot_triplet)) + rot_triplet=np.array(rot_triplet,dtype=int) + if (rot_triplet.size==0): + nRotTriplets=0 + nb=0 + #return np.array([]),0,0; + else: + nRotTriplets,nb = rot_triplet.shape; + + #print(nRotTriplets,nb,rot_triplet.flatten()) + + if (nb != 3 and nRotTriplets != 0 and ntot!= 0): + print('**ERROR: the number of column vectors in the rotating triplet must equal 3, the num of blades'); + new_seq = np.range(1,ntot) + else: + non_rotating = np.ones(ntot,dtype=int); + #print(non_rotating) + non_rotating[rot_triplet.flatten()] = 0; # if they are rotating, set them false; + a=np.array(np.nonzero(non_rotating)).flatten() + b=(rot_triplet.reshape(nRotTriplets*nb, 1)).flatten() + new_seq = np.concatenate((a,b)); + + #print(new_seq) + return new_seq,nRotTriplets,nb + +def eiganalysis(A, ndof2, ndof1): + # this line has to be the first line for this function + # to get the number of function arguments + nargin=len(locals()) + + mbc={} + m, ns = A.shape; + if(m!=ns): + print('**ERROR: the state-space matrix is not a square matrix.'); + + if nargin == 1: + ndof1 = 0; + ndof2 = ns/2; + + if np.mod(ns,2) != 0: + print('**ERROR: the input matrix is not of even order.'); + elif nargin == 2: + ndof1 = ns - 2*ndof2; + if ndof1 < 0: + print('**ERROR: ndof2 must be no larger than half the dimension of the state-space matrix.'); + else: + if ns != 2*ndof2 + ndof1: + print('**ERROR: the dimension of the state-space matrix must equal 2*ndof2 + ndof1.'); + + ndof = ndof2 + ndof1; + + origEvals, origEigenVects = np.linalg.eig(A); #,'nobalance' + # errorInSolution = norm(A * mbc.EigenVects - mbc.EigenVects* diag(mbc.EigenVals) ) + # these eigenvalues aren't sorted, so we just take the ones with + # positive imaginary parts to get the pairs for modes with damping < 1: + positiveImagEvals = np.argwhere( origEvals.imag > 0.0); + + mbc['Evals'] = origEvals[positiveImagEvals]; + row=np.array(list(range(0,ndof2))+list(range(ndof2*2+1,ns))) + col=positiveImagEvals + + mbc['EigenVects']=origEigenVects[row,col].transpose(); # save q2 and q1, throw away q2_dot + # print('GGGGGG ', row.shape,col.shape) + # with open('AAAA', "a") as f: + # np.savetxt(f,mbc['EigenVects'].imag,fmt='%5.4f') + # f.write('\n') + + EigenVects_save=origEigenVects[:,positiveImagEvals]; # save these for VTK visualization; + + #print('EigenVects_save shape ',EigenVects_save[:,:,0].shape, origEigenVects.shape) + + real_Evals = mbc['Evals'].real; + imag_Evals = mbc['Evals'].imag; + + mbc['NaturalFrequencies'] = np.sqrt( real_Evals**2 + imag_Evals**2 ); + mbc['DampRatios'] = -real_Evals/mbc['NaturalFrequencies'] + mbc['DampedFrequencies'] = imag_Evals; + + mbc['NumRigidBodyModes'] = ndof - len(positiveImagEvals); + + mbc['NaturalFreqs_Hz'] = mbc['NaturalFrequencies']/(2.0*np.pi) + mbc['DampedFreqs_Hz'] = mbc['DampedFrequencies']/(2.0*np.pi); + mbc['MagnitudeModes'] = np.abs(mbc['EigenVects']); + mbc['PhaseModes_deg'] = np.angle(mbc['EigenVects'])*180.0/np.pi; + + # print(real_Evals[0:10]) + # print(imag_Evals[0:10]) + # print(mbc['NaturalFrequencies'][0:10]) + # print(mbc['DampRatios'][0:10]) + # print(mbc['NaturalFreqs_Hz'][0:10]) + # print(mbc['DampedFreqs_Hz'][0:10]) + # print(mbc['MagnitudeModes'].shape) + # print(mbc['PhaseModes_deg'][0,0]) + + # with open('MagModes.txt', "a") as f: + # np.savetxt(f,mbc['MagnitudeModes'],fmt='%g') + # f.write('\n') + + # with open('PhsModes.txt', "a") as f: + # np.savetxt(f,mbc['PhaseModes_deg'],fmt='%g') + # f.write('\n') + + # with open('NFreqs_Hz.txt', "a") as f: + # np.savetxt(f,mbc['NaturalFreqs_Hz'],fmt='%g') + # f.write('\n') + # with open('DampRatios.txt', "a") as f: + # np.savetxt(f,mbc['DampRatios'],fmt='%g') + # f.write('\n') + # with open('DFreqs_Hz.txt', "a") as f: + # np.savetxt(f,mbc['DampedFreqs_Hz'],fmt='%g') + # f.write('\n') + + return mbc,EigenVects_save[:,:,0] + +def fx_mbc3(FileNames, verbose=True, removeTwrAzimuth=False): + """ + Original contribution by: Srinivasa B. Ramisett, ramisettisrinivas@yahoo.com, http://ramisetti.github.io + """ + MBC={} + matData, FAST_linData = gm.get_Mats(FileNames, verbose=verbose, removeTwrAzimuth=removeTwrAzimuth) + + # print('matData[Omega] ', matData['Omega']) + # print('matData[OmegaDot] ', matData['OmegaDot']) + + MBC['DescStates'] = matData['DescStates'] # save this in the MBC type for possible campbell_diagram processing later + MBC['ndof2'] = matData['ndof2'] + MBC['ndof1'] = matData['ndof1'] + MBC['RotSpeed_rpm'] = np.mean(matData['Omega'])*(30/np.pi); #rad/s to rpm + MBC['WindSpeed'] = np.mean(matData['WindSpeed']) # NOTE: might be NaN for old files + + # print('RotSpeed_rpm ',MBC['RotSpeed_rpm']) + # print('ndof1 ', MBC['ndof1']) + # print('ndof2 ', MBC['ndof2']) + # print(matData['RotTripletIndicesStates2']) + + # nb = 3; % number of blades required for MBC3 + # ---------- Multi-Blade-Coordinate transformation ------------------------------------------- + new_seq_dof2, dummy, nb = get_new_seq(matData['RotTripletIndicesStates2'],matData['ndof2']); # these are the first ndof2 states (not "first time derivative" states); these values are used to calculate matrix transformations + new_seq_dof1, dummy, nb2 = get_new_seq(matData['RotTripletIndicesStates1'],matData['ndof1']); # these are the first-order ndof1 states; these values are used to calculate matrix transformations + + # print('new_seq_dof2 ', new_seq_dof2) + # print('new_seq_dof1 ', new_seq_dof1) + # print('dummy ', dummy, ' nb ', nb) + #nb = max(nb,nb2); + #if (nb==0): + # print('*** fx_mbc3: no states were found, so assuming turbine has 3 blades. ***') + # nb = 3 + + + new_seq_states=np.concatenate((new_seq_dof2, new_seq_dof2+matData['ndof2'])) + if new_seq_dof1.size!=0: + new_seq_states=np.concatenate((new_seq_states,new_seq_dof1+matData['NumStates2'])) + + #new_seq_states = [new_seq_dof2; new_seq_dof2+matData['ndof2']; new_seq_dof1+matData['NumStates2']]; # combine the second-order states, including "first time derivatives", with first-order states (assumes ordering of displacements and velocities in state matrices); these values are used to calculate matrix transformations + # second-order NonRotating q2, second-order Rotating q2, + # second-order NonRotating q2_dot, second-order Rotating q2_dot, + # first-order NonRotating q1, first-order Rotating q1 + + + if nb == 3: + MBC['performedTransformation'] = True; + + if matData['n_RotTripletStates2'] + matData['n_RotTripletStates1'] < 1: + print('*** There are no rotating states. MBC transformation, therefore, cannot be performed.'); + # perhaps just warn and perform eigenanalysis anyway? + elif (matData['n_RotTripletStates2']*nb > matData['ndof2']): + print('**ERROR: the rotating second-order dof exceeds the total num of second-order dof'); + elif (matData['n_RotTripletStates1']*nb > matData['ndof1']): + print('**ERROR: the rotating first-order dof exceeds the total num of first-order dof'); + + new_seq_inp,dummy,dummy = get_new_seq(matData['RotTripletIndicesCntrlInpt'],matData['NumInputs']); + new_seq_out,dummy,dummy = get_new_seq(matData['RotTripletIndicesOutput'],matData['NumOutputs']); + + n_FixFrameStates2 = matData['ndof2'] - matData['n_RotTripletStates2']*nb; # fixed-frame second-order dof + n_FixFrameStates1 = matData['ndof1'] - matData['n_RotTripletStates1']*nb; # fixed-frame first-order dof + n_FixFrameInputs = matData['NumInputs'] - matData['n_RotTripletInputs']*nb; # fixed-frame control inputs + n_FixFrameOutputs = matData['NumOutputs'] - matData['n_RotTripletOutputs']*nb; # fixed-frame outputs + + #print(n_FixFrameOutputs,n_FixFrameInputs, n_FixFrameStates1, n_FixFrameStates2) + + if ( len(matData['Omega']) != matData['NAzimStep']): + print('**ERROR: the size of Omega vector must equal matData.NAzimStep, the num of azimuth steps') + if ( len(matData['OmegaDot']) != matData['NAzimStep']): + print('**ERROR: the size of OmegaDot vector must equal matData.NAzimStep, the num of azimuth steps'); + + + nLin = matData['A'].shape[-1] + MBC['A'] = np.zeros(matData['A'].shape) + MBC['B'] = np.zeros((len(new_seq_states),len(new_seq_inp),matData['NAzimStep'])) + if 'C' in matData.keys(): + MBC['C']=np.zeros(matData['C'].shape) + else: + MBC['C']=np.zeros((0,0,nLin)) + if 'D' in matData.keys(): + MBC['D']=np.zeros(matData['D'].shape) + else: + MBC['D']=np.zeros((0,0,nLin)) + + # print('new_seq_inp ',new_seq_inp) + # print('new_seq_out ',new_seq_out) + # print('new_seq_states ', new_seq_states) + + # begin azimuth loop + for iaz in reversed(range(matData['NAzimStep'])): + #(loop backwards so we don't reallocate memory each time [i.e. variables with iaz index aren't getting larger each time]) + + temp=np.arange(nb) + # compute azimuth positions of blades: + az = matData['Azimuth'][iaz]*np.pi/180.0 + 2*np.pi/nb* temp ; # Eq. 1, azimuth in radians + + # get rotor speed squared + OmegaSquared = matData['Omega'][iaz]**2; + + #print(OmegaSquared) + + # compute transformation matrices + cos_col = np.cos(az); + sin_col = np.sin(az); + + tt=np.column_stack((np.ones(3),cos_col,sin_col)) # Eq. 9, t_tilde + ttv = get_tt_inverse(sin_col, cos_col); # inverse of tt (computed analytically in function below) + tt2 = np.column_stack((np.zeros(3), -sin_col, cos_col)) # Eq. 16 a, t_tilde_2 + tt3 = np.column_stack((np.zeros(3), -cos_col, -sin_col)) # Eq. 16 b, t_tilde_3 + + #--- + T1 = np.eye(n_FixFrameStates2); # Eq. 11 for second-order states only + #print('B ',T1, n_FixFrameStates2, matData['n_RotTripletStates2']) + for ii in range(matData['n_RotTripletStates2']): + T1 = scp.block_diag(T1,tt) + + T1v = np.eye(n_FixFrameStates2); # inverse of T1 + for ii in range(matData['n_RotTripletStates2']): + T1v = scp.block_diag(T1v, ttv); + + T2 = np.zeros([n_FixFrameStates2,n_FixFrameStates2]); # Eq. 14 for second-order states only + for ii in range(matData['n_RotTripletStates2']): + T2 = scp.block_diag(T2, tt2); + + #print('T1, T1v, T2 ',T1.shape, T1v.shape, T2.shape) + #--- + T1q = np.eye(n_FixFrameStates1); # Eq. 11 for first-order states (eq. 8 in MBC3 Update document) + for ii in range(matData['n_RotTripletStates1']): + T1q = scp.block_diag(T1q, tt); + + T1qv = np.eye(n_FixFrameStates1); # inverse of T1q + for ii in range(matData['n_RotTripletStates1']): + T1qv = scp.block_diag(T1qv, ttv); + + T2q = np.zeros([n_FixFrameStates1,n_FixFrameStates1]); # Eq. 14 for first-order states (eq. 9 in MBC3 Update document) + for ii in range(matData['n_RotTripletStates1']): + T2q = scp.block_diag(T2q, tt2); + + #print('T1q, T1qv, T2q ',T1q.shape, T1qv.shape, T2q.shape) + # T1qc = np.eye(matData.NumHDInputs); # inverse of T1q + + #--- + T3 = np.zeros([n_FixFrameStates2,n_FixFrameStates2]); # Eq. 15 + for ii in range(matData['n_RotTripletStates2']): + T3 = scp.block_diag(T3, tt3); + + #--- + T1c = np.eye(n_FixFrameInputs); # Eq. 21 + for ii in range(matData['n_RotTripletInputs']): + T1c = scp.block_diag(T1c, tt) + + T1ov = np.eye(n_FixFrameOutputs); # inverse of Tlo (Eq. 23) + for ii in range(matData['n_RotTripletOutputs']): + T1ov = scp.block_diag(T1ov, ttv); + + #print('T3, T1c, T1ov ',T3.shape, T1c.shape, T1ov.shape, matData['A'].shape) + # mbc transformation of first-order matrices + # if ( MBC.EqnsOrder == 1 ) # activate later + + #print('Before ',T1c) + + if 'A' in matData: + # Eq. 29 + L1=np.concatenate((T1, np.zeros([matData['ndof2'],matData['ndof2']]), np.zeros([matData['ndof2'], matData['ndof1']])), axis=1) + L2=np.concatenate((matData['Omega'][iaz]*T2,T1,np.zeros([matData['ndof2'], matData['ndof1']])), axis=1) + L3=np.concatenate((np.zeros([matData['ndof1'], matData['ndof2']]),np.zeros([matData['ndof1'], matData['ndof2']]), T1q), axis=1) + L=np.matmul(matData['A'][new_seq_states[:,None],new_seq_states,iaz],np.concatenate((L1,L2,L3),axis=0)) + + R1=np.concatenate((matData['Omega'][iaz]*T2, np.zeros([matData['ndof2'],matData['ndof2']]), np.zeros([matData['ndof2'], matData['ndof1']])), axis=1) + R2=np.concatenate((OmegaSquared*T3 + matData['OmegaDot'][iaz]*T2, 2*matData['Omega'][iaz]*T2, np.zeros([matData['ndof2'], matData['ndof1']])),axis=1) + R3=np.concatenate((np.zeros([matData['ndof1'], matData['ndof2']]), np.zeros([matData['ndof1'], matData['ndof2']]), matData['Omega'][iaz]*T2q), axis=1) + + R=np.concatenate((R1,R2,R3),axis=0) + + MBC['A'][new_seq_states[:,None],new_seq_states,iaz]=np.matmul(scp.block_diag(T1v, T1v, T1qv),(L-R)) + + # ffname='AAA'+str(iaz)+'.txt' + # with open(ffname, "a") as f: + # np.savetxt(f,MBC['A'][:,:,iaz],fmt='%5.4f') + # f.write('\n') + + if 'B' in matData: + # Eq. 30 + MBC['B'][new_seq_states[:,None],new_seq_inp,iaz]=np.matmul(np.matmul(scp.block_diag(T1v, T1v, T1qv), matData['B'][new_seq_states[:,None],new_seq_inp,iaz]),T1c) + + # ffname='BBB'+str(iaz)+'.txt' + # with open(ffname, "a") as f: + # np.savetxt(f,MBC['B'][:,:,iaz],fmt='%5.4f') + # f.write('\n') + + if 'C' in matData: + # Eq. 31 + + L1=np.concatenate((T1, np.zeros([matData['ndof2'],matData['ndof2']]), np.zeros([matData['ndof2'], matData['ndof1']])),axis=1) + L2=np.concatenate((matData['Omega'][iaz]*T2, T1, np.zeros([matData['ndof2'], matData['ndof1']])), axis=1) + L3=np.concatenate((np.zeros([matData['ndof1'], matData['ndof2']]), np.zeros([matData['ndof1'], matData['ndof2']]), T1q), axis=1) + + MBC['C'][new_seq_out[:,None], new_seq_states,iaz]=np.matmul(np.matmul(T1ov,matData['C'][new_seq_out[:,None],new_seq_states,iaz]),np.concatenate((L1,L2,L3),axis=0)) + + # ffname='CCC'+str(iaz)+'.txt' + # with open(ffname, "a") as f: + # np.savetxt(f,MBC['C'][:,:,iaz],fmt='%5.4f') + # f.write('\n') + + if 'D' in matData: + # Eq. 32 + MBC['D'][new_seq_out[:,None],new_seq_inp,iaz] = np.matmul(np.matmul(T1ov,matData['D'][new_seq_out[:,None],new_seq_inp,iaz]), T1c) + + # ffname='DDD'+str(iaz)+'.txt' + # with open(ffname, "a") as f: + # np.savetxt(f,MBC['D'][:,:,iaz],fmt='%5.4f') + # f.write('\n') + + # end # end of azimuth loop + else: + print(' fx_mbc3 WARNING: Number of blades is ', str(nb), ' not 3. MBC transformation was not performed.') + MBC['performedTransformation'] = False; + + # initialize matrices + if 'A' in matData: + MBC['A'] = matData['A'] # initalize matrix + if 'B' in matData: + MBC['B'] = matData['B'] # initalize matrix + if 'C' in matData: + MBC['C'] = matData['C'] # initalize matrix + if 'D' in matData: + MBC['D'] = matData['D'] # initalize matrix + + # ------------- Eigensolution and Azimuth Averages ------------------------- + if 'A' in MBC: + MBC['AvgA'] = np.mean(MBC['A'],axis=2); # azimuth-average of azimuth-dependent MBC.A matrices + MBC['eigSol'], EigenVects_save = eiganalysis(MBC['AvgA'],matData['ndof2'], matData['ndof1']); + + # ffname='AAA_avg'+'.txt' + # with open(ffname, "a") as f: + # np.savetxt(f,MBC['AvgA'],fmt='%5.4f') + # f.write('\n') + + + # save eigenvectors (doing inverse of MBC3) for VTK visualization in FAST + # if nargout > 3 or nargin > 1: + # [VTK] = GetDataForVTK(MBC, matData, nb, EigenVects_save); + # if nargin > 1 + # WriteDataForVTK(VTK, ModeVizFileName) + + if 'B' in MBC: + MBC['AvgB'] = np.mean(MBC['B'],axis=2); # azimuth-average of azimuth-dependent MBC.B matrices + # ffname='BBB_avg'+'.txt' + # with open(ffname, "a") as f: + # np.savetxt(f,MBC['AvgB'],fmt='%5.4f') + # f.write('\n') + + if 'C' in MBC: + MBC['AvgC'] = np.mean(MBC['C'],axis=2); # azimuth-average of azimuth-dependent MBC.C matrices + # ffname='CCC_avg'+'.txt' + # with open(ffname, "a") as f: + # np.savetxt(f,MBC['AvgC'],fmt='%5.4f') + # f.write('\n') + + if 'D' in MBC: + MBC['AvgD'] = np.mean(MBC['D'],axis=2); # azimuth-average of azimuth-dependent MBC.D matrices + # ffname='DDD_avg'+'.txt' + # with open(ffname, "a") as f: + # np.savetxt(f,MBC['AvgD'],fmt='%5.4f') + # f.write('\n') + + return MBC, matData, FAST_linData + + +def runMBC(FileNames, NLinTimes=None, removeTwrAzimuth=False): + ''' + FileNames are .fst files + ''' + CampbellData={} + HubRad=None;TipRad=None; + BladeLen=None; TowerHt=None + dataFound=False; + for i in range(len(FileNames)): + basename=os.path.splitext(os.path.basename(FileNames[i]))[0] + dirname = os.path.dirname(FileNames[i]) + with open(FileNames[i]) as f: + datafile = f.readlines() + + if dataFound==False: + for line in datafile: + if 'EDFile' in line: + EDFile=line.split()[0].replace('"','') + with open(os.path.join(dirname,EDFile)) as edfile: + for edline in edfile: + if 'HubRad' in edline: + HubRad=float(edline.split()[0]) + elif 'TipRad' in edline: + TipRad=float(edline.split()[0]) + elif 'TowerHt' in edline: + TowerHt=float(edline.split()[0]) + + if((TipRad!=None and HubRad!=None) or TowerHt!=None): + BladeLen=HubRad-TipRad + dataFound=True + if(TowerHt==None or BladeLen==None): + print('TowerHt and BladeLen are not available!'); + sys.exit() + + if NLinTimes==None: + found=False + for line in datafile: + if found==False and 'NLinTimes' in line: + NLinTimes=int(line.split()[0]) + found=True + if NLinTimes<1: + print('NLinTimes should be greater than 0!') + sys.exit() + + linFileNames=[os.path.join(dirname,f'{basename}.{x:d}.lin') for x in range(1,NLinTimes+1)] + + print('Processing ', FileNames[i], ' file!') + MBC_data,getMatData,FAST_linData=fx_mbc3(linFileNames, removeTwrAzimuth=removeTwrAzimuth) + print('Multi-Blade Coordinate transformation completed!'); + print(' '); + CampbellData[i]=campbell_diagram_data(MBC_data,BladeLen,TowerHt) + + return CampbellData + + + +if __name__=='__main__': + pass + + # FileNames=['5MW_Land_ModeShapes-1.fst', '5MW_Land_ModeShapes-2.fst', '5MW_Land_ModeShapes-3.fst', '5MW_Land_ModeShapes-6.fst', '5MW_Land_ModeShapes-7.fst']; + #FileNames=['5MW_Land_BD_Linear-1.fst', '5MW_Land_BD_Linear-2.fst', '5MW_Land_BD_Linear-3.fst', '5MW_Land_BD_Linear-6.fst', '5MW_Land_BD_Linear-7.fst']; + + #FileNames=['5MW_Land_BD_Linear-1.fst']; + + #FileNames=['DLC-1.1/5MW_Land_BD_Linear-7.1.lin', 'DLC-1.1/5MW_Land_BD_Linear-7.2.lin'] + #FileNames=['/Users/sramiset/Desktop/OpenFAST/5MW_Land_BD_Linear/5MW_Land_BD_Linear-1.1.lin','/Users/sramiset/Desktop/OpenFAST/5MW_Land_BD_Linear/5MW_Land_BD_Linear-1.2.lin'] + # CampbellData=runMBC(FileNames) + # print('Preparing campbell diagram data!'); + # # TO DO read x-axis for wind speed or rotor speed from csv file + # #op_csv=pd.read_csv('input.csv', sep=',') + # OP=[2,4,6,8,10] + + # modeID_table,modesDesc=IdentifyModes(CampbellData) + + # #print(modesDesc) + + # nModes=modeID_table.shape[0] + # nRuns=modeID_table.shape[1] + # cols=[item[0] for item in list(modesDesc.values())] + # #cols.append('1P');cols.append('3P');cols.append('6P') + # #cols.append('9P');cols.append('12P') + # frequency=pd.DataFrame(np.nan, index=np.arange(nRuns), columns=cols) + # dampratio=pd.DataFrame(np.nan, index=np.arange(nRuns), columns=cols) + # FreqPlotData=np.zeros((nRuns,nModes)) + # DampPlotData=np.zeros((nRuns,nModes)) + # for i in range(nRuns): + # for modeID in range(len(modesDesc)): # list of modes we want to identify + # idx=int(modeID_table[modeID,i]) + # FreqPlotData[i,modeID]=CampbellData[i]['Modes'][idx]['NaturalFreq_Hz'] + # DampPlotData[i,modeID]=CampbellData[i]['Modes'][idx]['DampingRatio'] + # #print(i,modeID,modesDesc[modeID][0],FreqPlotData[i,modeID]) + # frequency.iloc[i,:]=FreqPlotData[i,:] + # dampratio.iloc[i,:]=DampPlotData[i,:] + + # for i in range(len(OP)): + # # for 15 DOF + # frequency.index.values[i]=OP[i] + # dampratio.index.values[i]=OP[i] + + # # import openpyxl + # # xfile = openpyxl.load_workbook('/Users/sramiset/Desktop/OpenFAST/mbc3_py/CampbellDiagram_Template.xlsx') + + # pCD.plotCampbellData(OP,frequency,dampratio) + + # frequency['1P']=np.nan + # frequency['3P']=np.nan + # frequency['6P']=np.nan + # frequency['9P']=np.nan + # frequency['12P']=np.nan + + # print(nRuns) + # for i in range(nRuns): + # # for 1P,3P,6P,9P,and 12P harmonics + # tmp=OP[i]/60.0 + # print(i,tmp) + # LZ=15 + # frequency.iloc[i,LZ]=tmp + # frequency.iloc[i,LZ+1]=3*tmp + # frequency.iloc[i,LZ+2]=6*tmp + # frequency.iloc[i,LZ+3]=9*tmp + # frequency.iloc[i,LZ+4]=12*tmp + # print(frequency) + # frequency.transpose().to_excel(r'CampbellData.xlsx') diff --git a/rosco/toolbox/ofTools/case_gen/run_FAST.py b/rosco/toolbox/ofTools/case_gen/run_FAST.py index dda46de9..48503553 100644 --- a/rosco/toolbox/ofTools/case_gen/run_FAST.py +++ b/rosco/toolbox/ofTools/case_gen/run_FAST.py @@ -143,7 +143,7 @@ def run_FAST(self): # lib dir if not in_weis: # in ROSCO - dll_dir = os.path.join(self.rosco_dir, 'rosco/controller/build/') + dll_dir = os.path.join(self.rosco_dir, 'lib') else: dll_dir = os.path.join(self.rosco_dir,'../local/lib/') diff --git a/rosco/toolbox/turbine.py b/rosco/toolbox/turbine.py index b5d12799..73965285 100644 --- a/rosco/toolbox/turbine.py +++ b/rosco/toolbox/turbine.py @@ -643,7 +643,7 @@ def __init__(self,performance_table, pitch_initial_rad, TSR_initial): f_performance = interpolate.interp1d(TSR_initial,performance_beta_max,bounds_error='False',kind='quadratic') # interpolate function for Cx(tsr) values performance_fine = f_performance(TSR_fine_ind) # Cx values at fine pitch performance_max_ind = np.where(performance_fine == np.max(performance_fine)) # Find max performance at fine pitch - self.TSR_opt = float(TSR_fine[performance_max_ind[0]]) # TSR to maximize Cx at fine pitch + self.TSR_opt = float(TSR_fine[performance_max_ind[0]][0]) # TSR to maximize Cx at fine pitch def interp_surface(self,pitch,TSR): ''' diff --git a/rosco/toolbox/utilities.py b/rosco/toolbox/utilities.py index 30de9d6e..15e4d550 100644 --- a/rosco/toolbox/utilities.py +++ b/rosco/toolbox/utilities.py @@ -90,6 +90,7 @@ def write_DISCON(turbine, controller, param_file='DISCON.IN', txt_filename='Cp_C file.write('!------- SIMULATION CONTROL ------------------------------------------------------------\n') file.write('{0:<12d} ! LoggingLevel - {{0: write no debug files, 1: write standard output .dbg-file, 2: LoggingLevel 1 + ROSCO LocalVars (.dbg2) 3: LoggingLevel 2 + complete avrSWAP-array (.dbg3)}}\n'.format(int(rosco_vt['LoggingLevel']))) file.write('{} ! DT_Out - {{Time step to output .dbg* files, or 0 to match sampling period of OpenFAST}}\n'.format(rosco_vt['DT_Out'])) + file.write('{:<11d} ! Ext_Interface - ({})\n'.format(int(rosco_vt['Ext_Interface']), input_descriptions['Ext_Interface'])) file.write('{:<11d} ! Echo - ({})\n'.format(int(rosco_vt['Echo']), input_descriptions['Echo'])) file.write('\n') file.write('!------- CONTROLLER FLAGS -------------------------------------------------\n') diff --git a/setup.py b/setup.py index 8b83331a..bda4a7c9 100644 --- a/setup.py +++ b/setup.py @@ -9,50 +9,32 @@ # CONDITIONS OF ANY KIND, either express or implied. See the License for the # specific language governing permissions and limitations under the License. -"""The setup script.""" - -# This setup file was made to mimic the setup.py script for https://github.com/NREL/floris/ -# Accessed on January 9, 2020 - -# Note: To use the 'upload' functionality of this file, you must: -# $ pip install twine - -import io import os -import sys -from shutil import rmtree, copy -import glob import platform -from setuptools import find_packages, setup, Command, Extension +import shutil +import sysconfig +from setuptools import setup, Extension from setuptools.command.build_ext import build_ext -from setuptools.command.install import install as _install - -from io import open - -# Package meta-data. -NAME = 'rosco' -DESCRIPTION = 'A reference open source controller toolset for wind turbine applications.' -URL = 'https://github.com/NREL/ROSCO' -EMAIL = 'daniel.zalkind@nrel.gov' -AUTHOR = 'NREL, National Wind Technology Center' -REQUIRES_PYTHON = '>=3.8' -VERSION = '2.8.0' - -# These packages are required for all of the code to be executed. -# - Maybe you can get away with older versions... -REQUIRED = [ - 'matplotlib', - 'numpy', - #'pytest', - 'scipy', - 'pyYAML', - #'future', - 'pandas' -] +####### +# This forces wheels to be platform specific +from setuptools.dist import Distribution +from wheel.bdist_wheel import bdist_wheel as _bdist_wheel + +class bdist_wheel(_bdist_wheel): + def finalize_options(self): + _bdist_wheel.finalize_options(self) + self.root_is_pure = False + +class BinaryDistribution(Distribution): + """Distribution which always forces a binary package with platform name""" + def has_ext_modules(foo): + return True +####### # For the CMake Extensions this_directory = os.path.abspath(os.path.dirname(__file__)) +build_dir = os.path.join(this_directory, "build") class CMakeExtension(Extension): @@ -82,14 +64,21 @@ def build_extension(self, ext): raise RuntimeError('Cannot find CMake executable') # Refresh build directory - localdir = os.path.join(this_directory, 'rosco','controller','install') - os.makedirs(localdir, exist_ok=True) + #os.makedirs(localdir, exist_ok=True) cmake_args = ['-DBUILD_SHARED_LIBS=OFF'] cmake_args += ['-DCMAKE_Fortran_FLAGS=-ffree-line-length-0'] - cmake_args += ['-DCMAKE_INSTALL_PREFIX={}'.format(localdir)] + cmake_args += ['-DCMAKE_INSTALL_PREFIX={}'.format(this_directory)] + + # Help Cmake find libraries in python locations + python_root = os.path.dirname( os.path.dirname( sysconfig.get_path('stdlib') ) ) + user_root = sysconfig.get_config_var("userbase") + cmake_args += [f'-DCMAKE_PREFIX_PATH={python_root}'] if platform.system() == 'Windows': + if not "FC" in os.environ: + os.environ["FC"] = "gfortran" + if "gfortran" in os.environ["FC"].lower(): cmake_args += ['-G', 'MinGW Makefiles'] elif self.compiler.compiler_type == 'msvc': @@ -97,100 +86,21 @@ def build_extension(self, ext): else: raise ValueError("Unable to find the system's Fortran compiler.") - self.build_temp = os.path.join( os.path.dirname( os.path.realpath(__file__) ), 'rosco', 'controller', 'build') - os.makedirs(localdir, exist_ok=True) + self.build_temp = build_dir + # Need fresh build directory for CMake os.makedirs(self.build_temp, exist_ok=True) - self.spawn(['cmake', '-S', ext.sourcedir, '-B', self.build_temp] + cmake_args) + self.spawn(['cmake', '-B', self.build_temp, '-S', ext.sourcedir] + cmake_args) self.spawn(['cmake', '--build', self.build_temp, '--target', 'install', '--config', 'Release']) - - -# All of the extensions -roscoExt = CMakeExtension('rosco',os.path.join('rosco','controller')) - -# The rest you shouldn't have to touch too much :) -# ------------------------------------------------ -# Except, perhaps the License and Trove Classifiers! -# If you do change the License, remember to change the Trove Classifier for that! - -here = os.path.abspath(os.path.dirname(__file__)) - -# Import the README and use it as the long-description. -try: - with io.open(os.path.join(here, 'README.md'), encoding='utf-8') as f: - long_description = '\n' + f.read() -except FileNotFoundError: - long_description = DESCRIPTION - -# Load the package's __version__.py module as a dictionary. -about = {} -if not VERSION: - project_slug = NAME.lower().replace("-", "_").replace(" ", "_") - with open(os.path.join(here, project_slug, '__version__.py')) as f: - exec(f.read(), about) -else: - about['__version__'] = VERSION - - -class UploadCommand(Command): - """Support setup.py upload.""" - - description = 'Build and publish the package.' - user_options = [] - - @staticmethod - def status(s): - """Prints things in bold.""" - print('\033[1m{0}\033[0m'.format(s)) - - def initialize_options(self): - pass - - def finalize_options(self): - pass - - def run(self): - try: - self.status('Removing previous builds...') - rmtree(os.path.join(here, 'dist')) - except OSError: - pass - - self.status('Building Source and Wheel (universal) distribution...') - os.system('{0} setup.py sdist bdist_wheel --universal'.format(sys.executable)) - - self.status('Uploading the package to PyPI via Twine...') - os.system('twine upload dist/*') - - self.status('Pushing git tags...') - os.system('git tag v{0}'.format(about['__version__'])) - os.system('git push --tags') - - sys.exit() - - - -metadata = dict( - name = NAME, - version = about['__version__'], - description = DESCRIPTION, - long_description = long_description, - long_description_content_type = 'text/markdown', - author = AUTHOR, - author_email = EMAIL, - url = URL, - install_requires = REQUIRED, - python_requires = REQUIRES_PYTHON, - packages = find_packages(exclude=["tests", "*.tests", "*.tests.*", "tests.*"]), - package_data = {'': ['*.yaml']}, - license = 'Apache License, Version 2.0', - cmdclass = {'build_ext': CMakeBuildExt, 'upload': UploadCommand}, - zip_safe = False, -) -if "--compile-rosco" in sys.argv: - metadata['ext_modules'] = [roscoExt] - sys.argv.remove("--compile-rosco") - -setup(**metadata) + +if __name__ == "__main__": + # Start with clean build directory + shutil.rmtree(build_dir, ignore_errors=True) + + setup(cmdclass={'bdist_wheel': bdist_wheel, 'build_ext': CMakeBuildExt}, + distclass=BinaryDistribution, + ext_modules=[ CMakeExtension('rosco',os.path.join('rosco','controller')) ], + ) +