diff --git a/.coveragerc b/.coveragerc index b602a0c6..16c960e0 100644 --- a/.coveragerc +++ b/.coveragerc @@ -1,5 +1,7 @@ [run] -omit = +omit = ard/utils/test_utils.py ard/farm_aero/placeholder.py - + flowers/visualization.py + flowers/optimization/optimization_interface.py + flowers/optimization/model_interface.py \ No newline at end of file diff --git a/.github/workflows/python-tests-consolidated.yaml b/.github/workflows/python-tests-consolidated.yaml index fb420d8c..5ff19dc7 100644 --- a/.github/workflows/python-tests-consolidated.yaml +++ b/.github/workflows/python-tests-consolidated.yaml @@ -46,8 +46,8 @@ jobs: run: | pip install .[dev] - test-unit: - name: Run unit tests + test-ard-unit: + name: Run Ard unit tests needs: setup-install strategy: matrix: @@ -86,9 +86,9 @@ jobs: run: | pytest --cov=ard --cov-fail-under=80 test/ard/unit - test-system: - name: Run system tests - needs: test-unit + test-ard-system: + name: Run Ard system tests + needs: test-ard-unit strategy: matrix: python-version: [3.12] # ["3.10", "3.11", "3.12", "3.13"] @@ -127,9 +127,90 @@ jobs: pytest --cov=ard --cov-fail-under=50 test/ard/system # pytest --cov=ard --cov-fail-under=80 test/ard/system + test-flowers-unit: + name: Run FLOWERS unit tests + needs: setup-install + strategy: + matrix: + python-version: [3.12] # ["3.10", "3.11", "3.12", "3.13"] + os: [macos-latest, ubuntu-latest, windows-latest] + include: + - os: ubuntu-latest + path: ~/.cache/pip + - os: macos-latest + path: ~/Library/Caches/pip + - os: windows-latest + path: ~\AppData\Local\pip\Cache + runs-on: ${{ matrix.os }} + steps: + - name: Checkout code + uses: actions/checkout@v4 + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python-version }} + - name: Get cache dependencies + id: cache + uses: actions/cache@v4 + with: + path: ${{ matrix.path }} + key: ${{ runner.os }}-pip-${{ hashFiles('requirements.txt') }} + restore-keys: + ${{ runner.os }}-pip- + - name: Install dependencies + run: | + python -m pip install --upgrade pip + - name: Install Ard + run: | + pip install .[dev] + - name: Run unit tests with coverage + run: | + pytest --cov=flowers --cov-fail-under=80 test/flowers/unit + + test-flowers-system: + name: Run FLOWERS system tests + needs: test-flowers-unit + strategy: + matrix: + python-version: [3.12] # ["3.10", "3.11", "3.12", "3.13"] + os: [macos-latest, ubuntu-latest, windows-latest] + include: + - os: ubuntu-latest + path: ~/.cache/pip + - os: macos-latest + path: ~/Library/Caches/pip + - os: windows-latest + path: ~\AppData\Local\pip\Cache + runs-on: ${{ matrix.os }} + steps: + - name: Checkout code + uses: actions/checkout@v4 + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python-version }} + - name: Get cache dependencies + id: cache + uses: actions/cache@v4 + with: + path: ${{ matrix.path }} + key: ${{ runner.os }}-pip-${{ hashFiles('requirements.txt') }} + restore-keys: + ${{ runner.os }}-pip- + - name: Install dependencies + run: | + python -m pip install --upgrade pip + - name: Install Ard + run: | + pip install .[dev] + - name: Run system tests with coverage + run: | + pytest --cov=flowers --cov-fail-under=50 test/flowers/system + # pytest --cov=flowers --cov-fail-under=80 test/flowers/system + find-examples: name: Find all examples - needs: [test-unit, test-system] + needs: [test-ard-unit, test-ard-system, test-flowers-unit, test-flowers-system] runs-on: ubuntu-latest steps: - name: Checkout code @@ -142,7 +223,7 @@ jobs: id: find_scripts run: | echo "scripts<> $GITHUB_OUTPUT - find examples -mindepth 2 -maxdepth 2 -name "*.py" | jq -R -s -c 'split("\n")[:-1]' >> $GITHUB_OUTPUT + find examples -mindepth 3 -maxdepth 3 -name "*.py" | jq -R -s -c 'split("\n")[:-1]' >> $GITHUB_OUTPUT echo "EOF" >> $GITHUB_OUTPUT outputs: scripts: ${{ steps.find_scripts.outputs.scripts }} @@ -207,7 +288,7 @@ jobs: find-examples_nb: name: Find all example notebooks - needs: [test-unit, test-system] + needs: [test-ard-unit, test-ard-system, test-flowers-unit, test-flowers-system] if: needs.find-examples_nb.outputs.scripts != '[]' # skip empty runs-on: ubuntu-latest steps: @@ -221,10 +302,10 @@ jobs: id: find_notebooks run: | echo "scripts<> $GITHUB_OUTPUT - find examples -mindepth 2 -maxdepth 2 -name "*.ipynb" | jq -R -s -c 'split("\n")[:-1]' >> $GITHUB_OUTPUT + find examples -mindepth 3 -maxdepth 3 -name "*.ipynb" | jq -R -s -c 'split("\n")[:-1]' >> $GITHUB_OUTPUT echo "EOF" >> $GITHUB_OUTPUT outputs: scripts: ${{ steps.find_notebooks.outputs.scripts }} diff --git a/.gitignore b/.gitignore index 20e9a213..3870c2c7 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,5 @@ -### ARD DEVELOPMENT IGNORES +### ARD DEVELOPMENT IGNORES .vscode case_files diff --git a/assets/logomaker/inputs/windio.yaml b/assets/logomaker/inputs/windio.yaml index 83e0cdb0..e076be17 100644 --- a/assets/logomaker/inputs/windio.yaml +++ b/assets/logomaker/inputs/windio.yaml @@ -3303,7 +3303,7 @@ site: - 511.7766126774107 energy_resource: name: Example energy resource - wind_resource: !include ../../../examples/data/windIO-plant_wind-resource_wrg-example.yaml + wind_resource: !include ../../../examples/ard/data/windIO-plant_wind-resource_wrg-example.yaml wind_farm: name: LogoFarm layouts: @@ -3322,7 +3322,7 @@ wind_farm: 1250.0, 1250.0, 1250.0, 1250.0, 1250.0, 2500.0, 2500.0, 2500.0, 2500.0, 2500.0 ] - turbine: !include ../../../examples/data/windIO-plant_turbine_IEA-3.4MW-130m-RWT.yaml + turbine: !include ../../../examples/ard/data/windIO-plant_turbine_IEA-3.4MW-130m-RWT.yaml electrical_substations: - electrical_substation: coordinates: diff --git a/examples/01_onshore/inputs/ard_system.yaml b/examples/ard/01_onshore/inputs/ard_system.yaml similarity index 100% rename from examples/01_onshore/inputs/ard_system.yaml rename to examples/ard/01_onshore/inputs/ard_system.yaml diff --git a/examples/01_onshore/inputs/windio.yaml b/examples/ard/01_onshore/inputs/windio.yaml similarity index 100% rename from examples/01_onshore/inputs/windio.yaml rename to examples/ard/01_onshore/inputs/windio.yaml diff --git a/examples/01_onshore/optimization_demo.ipynb b/examples/ard/01_onshore/optimization_demo.ipynb similarity index 100% rename from examples/01_onshore/optimization_demo.ipynb rename to examples/ard/01_onshore/optimization_demo.ipynb diff --git a/examples/02_offshore_fixed/inputs/ard_system.yaml b/examples/ard/02_offshore_fixed/inputs/ard_system.yaml similarity index 100% rename from examples/02_offshore_fixed/inputs/ard_system.yaml rename to examples/ard/02_offshore_fixed/inputs/ard_system.yaml diff --git a/examples/02_offshore_fixed/inputs/windio.yaml b/examples/ard/02_offshore_fixed/inputs/windio.yaml similarity index 100% rename from examples/02_offshore_fixed/inputs/windio.yaml rename to examples/ard/02_offshore_fixed/inputs/windio.yaml diff --git a/examples/02_offshore_fixed/optimization_demo.ipynb b/examples/ard/02_offshore_fixed/optimization_demo.ipynb similarity index 100% rename from examples/02_offshore_fixed/optimization_demo.ipynb rename to examples/ard/02_offshore_fixed/optimization_demo.ipynb diff --git a/examples/03_offshore_floating_custom_system/inputs/ard_system.yaml b/examples/ard/03_offshore_floating_custom_system/inputs/ard_system.yaml similarity index 100% rename from examples/03_offshore_floating_custom_system/inputs/ard_system.yaml rename to examples/ard/03_offshore_floating_custom_system/inputs/ard_system.yaml diff --git a/examples/03_offshore_floating_custom_system/inputs/windio.yaml b/examples/ard/03_offshore_floating_custom_system/inputs/windio.yaml similarity index 100% rename from examples/03_offshore_floating_custom_system/inputs/windio.yaml rename to examples/ard/03_offshore_floating_custom_system/inputs/windio.yaml diff --git a/examples/03_offshore_floating_custom_system/optimization_demo.ipynb b/examples/ard/03_offshore_floating_custom_system/optimization_demo.ipynb similarity index 100% rename from examples/03_offshore_floating_custom_system/optimization_demo.ipynb rename to examples/ard/03_offshore_floating_custom_system/optimization_demo.ipynb diff --git a/examples/05_onshore_batch/inputs/NREL_Reference_5MW_126.csv b/examples/ard/05_onshore_batch/inputs/NREL_Reference_5MW_126.csv similarity index 100% rename from examples/05_onshore_batch/inputs/NREL_Reference_5MW_126.csv rename to examples/ard/05_onshore_batch/inputs/NREL_Reference_5MW_126.csv diff --git a/examples/05_onshore_batch/inputs/ard_system.yaml b/examples/ard/05_onshore_batch/inputs/ard_system.yaml similarity index 100% rename from examples/05_onshore_batch/inputs/ard_system.yaml rename to examples/ard/05_onshore_batch/inputs/ard_system.yaml diff --git a/examples/05_onshore_batch/inputs/h2i/driver_config.yaml b/examples/ard/05_onshore_batch/inputs/h2i/driver_config.yaml similarity index 100% rename from examples/05_onshore_batch/inputs/h2i/driver_config.yaml rename to examples/ard/05_onshore_batch/inputs/h2i/driver_config.yaml diff --git a/examples/05_onshore_batch/inputs/h2i/h2i_ard.yaml b/examples/ard/05_onshore_batch/inputs/h2i/h2i_ard.yaml similarity index 100% rename from examples/05_onshore_batch/inputs/h2i/h2i_ard.yaml rename to examples/ard/05_onshore_batch/inputs/h2i/h2i_ard.yaml diff --git a/examples/05_onshore_batch/inputs/h2i/plant_config.yaml b/examples/ard/05_onshore_batch/inputs/h2i/plant_config.yaml similarity index 100% rename from examples/05_onshore_batch/inputs/h2i/plant_config.yaml rename to examples/ard/05_onshore_batch/inputs/h2i/plant_config.yaml diff --git a/examples/05_onshore_batch/inputs/open-meteo-56.20N8.54E86m-short.yaml b/examples/ard/05_onshore_batch/inputs/open-meteo-56.20N8.54E86m-short.yaml similarity index 100% rename from examples/05_onshore_batch/inputs/open-meteo-56.20N8.54E86m-short.yaml rename to examples/ard/05_onshore_batch/inputs/open-meteo-56.20N8.54E86m-short.yaml diff --git a/examples/05_onshore_batch/inputs/open-meteo-56.20N8.54E86m.yaml b/examples/ard/05_onshore_batch/inputs/open-meteo-56.20N8.54E86m.yaml similarity index 100% rename from examples/05_onshore_batch/inputs/open-meteo-56.20N8.54E86m.yaml rename to examples/ard/05_onshore_batch/inputs/open-meteo-56.20N8.54E86m.yaml diff --git a/examples/05_onshore_batch/inputs/power_thrust_table_ccblade_NREL-5p0-126-RWT.csv b/examples/ard/05_onshore_batch/inputs/power_thrust_table_ccblade_NREL-5p0-126-RWT.csv similarity index 100% rename from examples/05_onshore_batch/inputs/power_thrust_table_ccblade_NREL-5p0-126-RWT.csv rename to examples/ard/05_onshore_batch/inputs/power_thrust_table_ccblade_NREL-5p0-126-RWT.csv diff --git a/examples/05_onshore_batch/inputs/resource_prep/convert_weather_to_speed_dir.py b/examples/ard/05_onshore_batch/inputs/resource_prep/convert_weather_to_speed_dir.py similarity index 100% rename from examples/05_onshore_batch/inputs/resource_prep/convert_weather_to_speed_dir.py rename to examples/ard/05_onshore_batch/inputs/resource_prep/convert_weather_to_speed_dir.py diff --git a/examples/05_onshore_batch/inputs/resource_prep/open-meteo-56.20N8.54E86m.csv b/examples/ard/05_onshore_batch/inputs/resource_prep/open-meteo-56.20N8.54E86m.csv similarity index 100% rename from examples/05_onshore_batch/inputs/resource_prep/open-meteo-56.20N8.54E86m.csv rename to examples/ard/05_onshore_batch/inputs/resource_prep/open-meteo-56.20N8.54E86m.csv diff --git a/examples/05_onshore_batch/inputs/windIO-plant_turbine_NREL-5.0MW-126m-RWT.yaml b/examples/ard/05_onshore_batch/inputs/windIO-plant_turbine_NREL-5.0MW-126m-RWT.yaml similarity index 100% rename from examples/05_onshore_batch/inputs/windIO-plant_turbine_NREL-5.0MW-126m-RWT.yaml rename to examples/ard/05_onshore_batch/inputs/windIO-plant_turbine_NREL-5.0MW-126m-RWT.yaml diff --git a/examples/05_onshore_batch/inputs/windio.yaml b/examples/ard/05_onshore_batch/inputs/windio.yaml similarity index 100% rename from examples/05_onshore_batch/inputs/windio.yaml rename to examples/ard/05_onshore_batch/inputs/windio.yaml diff --git a/examples/05_onshore_batch/onshore-batch.ipynb b/examples/ard/05_onshore_batch/onshore-batch.ipynb similarity index 100% rename from examples/05_onshore_batch/onshore-batch.ipynb rename to examples/ard/05_onshore_batch/onshore-batch.ipynb diff --git a/examples/05_onshore_batch/run_wind_ard.py b/examples/ard/05_onshore_batch/run_wind_ard.py similarity index 100% rename from examples/05_onshore_batch/run_wind_ard.py rename to examples/ard/05_onshore_batch/run_wind_ard.py diff --git a/examples/06_onshore_multiobjective/inputs/ard_system.yaml b/examples/ard/06_onshore_multiobjective/inputs/ard_system.yaml similarity index 100% rename from examples/06_onshore_multiobjective/inputs/ard_system.yaml rename to examples/ard/06_onshore_multiobjective/inputs/ard_system.yaml diff --git a/examples/06_onshore_multiobjective/inputs/windio.yaml b/examples/ard/06_onshore_multiobjective/inputs/windio.yaml similarity index 100% rename from examples/06_onshore_multiobjective/inputs/windio.yaml rename to examples/ard/06_onshore_multiobjective/inputs/windio.yaml diff --git a/examples/06_onshore_multiobjective/optimization_demo.ipynb b/examples/ard/06_onshore_multiobjective/optimization_demo.ipynb similarity index 100% rename from examples/06_onshore_multiobjective/optimization_demo.ipynb rename to examples/ard/06_onshore_multiobjective/optimization_demo.ipynb diff --git a/examples/data/offshore/GulfOfMaine_bathymetry_100x99.txt b/examples/ard/data/offshore/GulfOfMaine_bathymetry_100x99.txt similarity index 100% rename from examples/data/offshore/GulfOfMaine_bathymetry_100x99.txt rename to examples/ard/data/offshore/GulfOfMaine_bathymetry_100x99.txt diff --git a/examples/data/windIO-plant_turbine_IEA-22MW-284m-RWT.yaml b/examples/ard/data/windIO-plant_turbine_IEA-22MW-284m-RWT.yaml similarity index 100% rename from examples/data/windIO-plant_turbine_IEA-22MW-284m-RWT.yaml rename to examples/ard/data/windIO-plant_turbine_IEA-22MW-284m-RWT.yaml diff --git a/examples/data/windIO-plant_turbine_IEA-3.4MW-130m-RWT.yaml b/examples/ard/data/windIO-plant_turbine_IEA-3.4MW-130m-RWT.yaml similarity index 100% rename from examples/data/windIO-plant_turbine_IEA-3.4MW-130m-RWT.yaml rename to examples/ard/data/windIO-plant_turbine_IEA-3.4MW-130m-RWT.yaml diff --git a/examples/data/windIO-plant_wind-resource_wrg-example.yaml b/examples/ard/data/windIO-plant_wind-resource_wrg-example.yaml similarity index 100% rename from examples/data/windIO-plant_wind-resource_wrg-example.yaml rename to examples/ard/data/windIO-plant_wind-resource_wrg-example.yaml diff --git a/examples/flowers/_archive/compute_AEP.py b/examples/flowers/_archive/compute_AEP.py new file mode 100644 index 00000000..836248cb --- /dev/null +++ b/examples/flowers/_archive/compute_AEP.py @@ -0,0 +1,23 @@ +import numpy as np +import pandas as pd + +from flowers import FlowersModel + +# Create generic layout +D = 126.0 +layout_x = D * np.array( + [0.0, 0.0, 0.0, 7.0, 7.0, 7.0, 14.0, 14.0, 14.0, 21.0, 21.0, 21.0, 28.0, 28.0, 28.0] +) +layout_y = D * np.array( + [0.0, 7.0, 14.0, 0.0, 7.0, 14.0, 0.0, 7.0, 14.0, 0.0, 7.0, 14.0, 0.0, 7.0, 14.0] +) + +# Load in wind data +df = pd.read_csv("../data/HKW_wind_rose.csv") + +# Setup FLOWERS model +flowers_model = FlowersModel(df, layout_x, layout_y) + +# Calculate the AEP +aep = flowers_model.calculate_aep() +print(aep) diff --git a/examples/flowers/data/HKW_wind_rose.csv b/examples/flowers/data/HKW_wind_rose.csv new file mode 100644 index 00000000..f608c6aa --- /dev/null +++ b/examples/flowers/data/HKW_wind_rose.csv @@ -0,0 +1,1873 @@ +wd,ws,freq_val +0.0000000000,0.0000000000,0.0000228154 +0.0000000000,1.0000000000,0.0002224504 +0.0000000000,2.0000000000,0.0003650468 +0.0000000000,3.0000000000,0.0004848277 +0.0000000000,4.0000000000,0.0008726899 +0.0000000000,5.0000000000,0.0008669861 +0.0000000000,6.0000000000,0.0009753593 +0.0000000000,7.0000000000,0.0009810632 +0.0000000000,8.0000000000,0.0009354324 +0.0000000000,9.0000000000,0.0009753593 +0.0000000000,10.0000000000,0.0007300935 +0.0000000000,11.0000000000,0.0005817933 +0.0000000000,12.0000000000,0.0005190509 +0.0000000000,13.0000000000,0.0003935661 +0.0000000000,14.0000000000,0.0003080082 +0.0000000000,15.0000000000,0.0002509697 +0.0000000000,16.0000000000,0.0001882272 +0.0000000000,17.0000000000,0.0000741501 +0.0000000000,18.0000000000,0.0000171116 +0.0000000000,19.0000000000,0.0000228154 +0.0000000000,20.0000000000,0.0000000000 +0.0000000000,21.0000000000,0.0000057039 +0.0000000000,22.0000000000,0.0000000000 +0.0000000000,23.0000000000,0.0000000000 +0.0000000000,24.0000000000,0.0000000000 +0.0000000000,25.0000000000,0.0000000000 +5.0000000000,0.0000000000,0.0000171116 +5.0000000000,1.0000000000,0.0001996350 +5.0000000000,2.0000000000,0.0004449008 +5.0000000000,3.0000000000,0.0006445357 +5.0000000000,4.0000000000,0.0007186858 +5.0000000000,5.0000000000,0.0007871321 +5.0000000000,6.0000000000,0.0008783938 +5.0000000000,7.0000000000,0.0011407712 +5.0000000000,8.0000000000,0.0011008442 +5.0000000000,9.0000000000,0.0010209902 +5.0000000000,10.0000000000,0.0007072781 +5.0000000000,11.0000000000,0.0006331280 +5.0000000000,12.0000000000,0.0004848277 +5.0000000000,13.0000000000,0.0003593429 +5.0000000000,14.0000000000,0.0003194159 +5.0000000000,15.0000000000,0.0002053388 +5.0000000000,16.0000000000,0.0001026694 +5.0000000000,17.0000000000,0.0000627424 +5.0000000000,18.0000000000,0.0000228154 +5.0000000000,19.0000000000,0.0000000000 +5.0000000000,20.0000000000,0.0000057039 +5.0000000000,21.0000000000,0.0000000000 +5.0000000000,22.0000000000,0.0000000000 +5.0000000000,23.0000000000,0.0000000000 +5.0000000000,24.0000000000,0.0000000000 +5.0000000000,25.0000000000,0.0000000000 +10.0000000000,0.0000000000,0.0000285193 +10.0000000000,1.0000000000,0.0002167465 +10.0000000000,2.0000000000,0.0003479352 +10.0000000000,3.0000000000,0.0005361624 +10.0000000000,4.0000000000,0.0006559434 +10.0000000000,5.0000000000,0.0007985398 +10.0000000000,6.0000000000,0.0008898015 +10.0000000000,7.0000000000,0.0009867671 +10.0000000000,8.0000000000,0.0011008442 +10.0000000000,9.0000000000,0.0009753593 +10.0000000000,10.0000000000,0.0007814282 +10.0000000000,11.0000000000,0.0006217203 +10.0000000000,12.0000000000,0.0005076432 +10.0000000000,13.0000000000,0.0003707506 +10.0000000000,14.0000000000,0.0002395619 +10.0000000000,15.0000000000,0.0001197810 +10.0000000000,16.0000000000,0.0000969655 +10.0000000000,17.0000000000,0.0000228154 +10.0000000000,18.0000000000,0.0000171116 +10.0000000000,19.0000000000,0.0000000000 +10.0000000000,20.0000000000,0.0000171116 +10.0000000000,21.0000000000,0.0000000000 +10.0000000000,22.0000000000,0.0000000000 +10.0000000000,23.0000000000,0.0000000000 +10.0000000000,24.0000000000,0.0000000000 +10.0000000000,25.0000000000,0.0000000000 +15.0000000000,0.0000000000,0.0000342231 +15.0000000000,1.0000000000,0.0001996350 +15.0000000000,2.0000000000,0.0004391969 +15.0000000000,3.0000000000,0.0004391969 +15.0000000000,4.0000000000,0.0006958704 +15.0000000000,5.0000000000,0.0008384668 +15.0000000000,6.0000000000,0.0009468401 +15.0000000000,7.0000000000,0.0010209902 +15.0000000000,8.0000000000,0.0010209902 +15.0000000000,9.0000000000,0.0008898015 +15.0000000000,10.0000000000,0.0007186858 +15.0000000000,11.0000000000,0.0006274241 +15.0000000000,12.0000000000,0.0004677162 +15.0000000000,13.0000000000,0.0002966005 +15.0000000000,14.0000000000,0.0001483003 +15.0000000000,15.0000000000,0.0000684463 +15.0000000000,16.0000000000,0.0000285193 +15.0000000000,17.0000000000,0.0000342231 +15.0000000000,18.0000000000,0.0000171116 +15.0000000000,19.0000000000,0.0000399270 +15.0000000000,20.0000000000,0.0000057039 +15.0000000000,21.0000000000,0.0000114077 +15.0000000000,22.0000000000,0.0000000000 +15.0000000000,23.0000000000,0.0000000000 +15.0000000000,24.0000000000,0.0000000000 +15.0000000000,25.0000000000,0.0000000000 +20.0000000000,0.0000000000,0.0000399270 +20.0000000000,1.0000000000,0.0001597080 +20.0000000000,2.0000000000,0.0003536391 +20.0000000000,3.0000000000,0.0004848277 +20.0000000000,4.0000000000,0.0006787588 +20.0000000000,5.0000000000,0.0007415013 +20.0000000000,6.0000000000,0.0011008442 +20.0000000000,7.0000000000,0.0010209902 +20.0000000000,8.0000000000,0.0011065480 +20.0000000000,9.0000000000,0.0010837326 +20.0000000000,10.0000000000,0.0007529090 +20.0000000000,11.0000000000,0.0005532740 +20.0000000000,12.0000000000,0.0004449008 +20.0000000000,13.0000000000,0.0003308236 +20.0000000000,14.0000000000,0.0001996350 +20.0000000000,15.0000000000,0.0001254848 +20.0000000000,16.0000000000,0.0000969655 +20.0000000000,17.0000000000,0.0000513347 +20.0000000000,18.0000000000,0.0000171116 +20.0000000000,19.0000000000,0.0000000000 +20.0000000000,20.0000000000,0.0000000000 +20.0000000000,21.0000000000,0.0000000000 +20.0000000000,22.0000000000,0.0000000000 +20.0000000000,23.0000000000,0.0000000000 +20.0000000000,24.0000000000,0.0000000000 +20.0000000000,25.0000000000,0.0000000000 +25.0000000000,0.0000000000,0.0000627424 +25.0000000000,1.0000000000,0.0002053388 +25.0000000000,2.0000000000,0.0003422313 +25.0000000000,3.0000000000,0.0005304586 +25.0000000000,4.0000000000,0.0007243897 +25.0000000000,5.0000000000,0.0008384668 +25.0000000000,6.0000000000,0.0008099475 +25.0000000000,7.0000000000,0.0009981748 +25.0000000000,8.0000000000,0.0011350673 +25.0000000000,9.0000000000,0.0009810632 +25.0000000000,10.0000000000,0.0008099475 +25.0000000000,11.0000000000,0.0006559434 +25.0000000000,12.0000000000,0.0004049738 +25.0000000000,13.0000000000,0.0002966005 +25.0000000000,14.0000000000,0.0002053388 +25.0000000000,15.0000000000,0.0001197810 +25.0000000000,16.0000000000,0.0000969655 +25.0000000000,17.0000000000,0.0000513347 +25.0000000000,18.0000000000,0.0000171116 +25.0000000000,19.0000000000,0.0000114077 +25.0000000000,20.0000000000,0.0000000000 +25.0000000000,21.0000000000,0.0000000000 +25.0000000000,22.0000000000,0.0000000000 +25.0000000000,23.0000000000,0.0000000000 +25.0000000000,24.0000000000,0.0000000000 +25.0000000000,25.0000000000,0.0000000000 +30.0000000000,0.0000000000,0.0000399270 +30.0000000000,1.0000000000,0.0001825234 +30.0000000000,2.0000000000,0.0003821583 +30.0000000000,3.0000000000,0.0004734200 +30.0000000000,4.0000000000,0.0005076432 +30.0000000000,5.0000000000,0.0007757244 +30.0000000000,6.0000000000,0.0009183208 +30.0000000000,7.0000000000,0.0010381018 +30.0000000000,8.0000000000,0.0009981748 +30.0000000000,9.0000000000,0.0010495095 +30.0000000000,10.0000000000,0.0008612822 +30.0000000000,11.0000000000,0.0007015743 +30.0000000000,12.0000000000,0.0005418663 +30.0000000000,13.0000000000,0.0003479352 +30.0000000000,14.0000000000,0.0002281542 +30.0000000000,15.0000000000,0.0001140771 +30.0000000000,16.0000000000,0.0000855578 +30.0000000000,17.0000000000,0.0000285193 +30.0000000000,18.0000000000,0.0000171116 +30.0000000000,19.0000000000,0.0000000000 +30.0000000000,20.0000000000,0.0000000000 +30.0000000000,21.0000000000,0.0000000000 +30.0000000000,22.0000000000,0.0000000000 +30.0000000000,23.0000000000,0.0000000000 +30.0000000000,24.0000000000,0.0000000000 +30.0000000000,25.0000000000,0.0000000000 +35.0000000000,0.0000000000,0.0000342231 +35.0000000000,1.0000000000,0.0001711157 +35.0000000000,2.0000000000,0.0003365275 +35.0000000000,3.0000000000,0.0005418663 +35.0000000000,4.0000000000,0.0005418663 +35.0000000000,5.0000000000,0.0006844627 +35.0000000000,6.0000000000,0.0009183208 +35.0000000000,7.0000000000,0.0009639516 +35.0000000000,8.0000000000,0.0009639516 +35.0000000000,9.0000000000,0.0009297285 +35.0000000000,10.0000000000,0.0008384668 +35.0000000000,11.0000000000,0.0007415013 +35.0000000000,12.0000000000,0.0006103126 +35.0000000000,13.0000000000,0.0004391969 +35.0000000000,14.0000000000,0.0002794889 +35.0000000000,15.0000000000,0.0002281542 +35.0000000000,16.0000000000,0.0000741501 +35.0000000000,17.0000000000,0.0000399270 +35.0000000000,18.0000000000,0.0000114077 +35.0000000000,19.0000000000,0.0000000000 +35.0000000000,20.0000000000,0.0000000000 +35.0000000000,21.0000000000,0.0000000000 +35.0000000000,22.0000000000,0.0000000000 +35.0000000000,23.0000000000,0.0000000000 +35.0000000000,24.0000000000,0.0000000000 +35.0000000000,25.0000000000,0.0000000000 +40.0000000000,0.0000000000,0.0000285193 +40.0000000000,1.0000000000,0.0001825234 +40.0000000000,2.0000000000,0.0003479352 +40.0000000000,3.0000000000,0.0004791239 +40.0000000000,4.0000000000,0.0003992699 +40.0000000000,5.0000000000,0.0007472051 +40.0000000000,6.0000000000,0.0008669861 +40.0000000000,7.0000000000,0.0010381018 +40.0000000000,8.0000000000,0.0010381018 +40.0000000000,9.0000000000,0.0011407712 +40.0000000000,10.0000000000,0.0009297285 +40.0000000000,11.0000000000,0.0008955054 +40.0000000000,12.0000000000,0.0008384668 +40.0000000000,13.0000000000,0.0006103126 +40.0000000000,14.0000000000,0.0003992699 +40.0000000000,15.0000000000,0.0002623774 +40.0000000000,16.0000000000,0.0001197810 +40.0000000000,17.0000000000,0.0000456308 +40.0000000000,18.0000000000,0.0000570386 +40.0000000000,19.0000000000,0.0000171116 +40.0000000000,20.0000000000,0.0000057039 +40.0000000000,21.0000000000,0.0000000000 +40.0000000000,22.0000000000,0.0000000000 +40.0000000000,23.0000000000,0.0000000000 +40.0000000000,24.0000000000,0.0000000000 +40.0000000000,25.0000000000,0.0000000000 +45.0000000000,0.0000000000,0.0000228154 +45.0000000000,1.0000000000,0.0001368925 +45.0000000000,2.0000000000,0.0003080082 +45.0000000000,3.0000000000,0.0004391969 +45.0000000000,4.0000000000,0.0006616473 +45.0000000000,5.0000000000,0.0007700205 +45.0000000000,6.0000000000,0.0008726899 +45.0000000000,7.0000000000,0.0010095825 +45.0000000000,8.0000000000,0.0011122519 +45.0000000000,9.0000000000,0.0013175907 +45.0000000000,10.0000000000,0.0010552133 +45.0000000000,11.0000000000,0.0010038786 +45.0000000000,12.0000000000,0.0008498745 +45.0000000000,13.0000000000,0.0006844627 +45.0000000000,14.0000000000,0.0004563085 +45.0000000000,15.0000000000,0.0002794889 +45.0000000000,16.0000000000,0.0001026694 +45.0000000000,17.0000000000,0.0000399270 +45.0000000000,18.0000000000,0.0000513347 +45.0000000000,19.0000000000,0.0000171116 +45.0000000000,20.0000000000,0.0000114077 +45.0000000000,21.0000000000,0.0000000000 +45.0000000000,22.0000000000,0.0000057039 +45.0000000000,23.0000000000,0.0000000000 +45.0000000000,24.0000000000,0.0000000000 +45.0000000000,25.0000000000,0.0000000000 +50.0000000000,0.0000000000,0.0000342231 +50.0000000000,1.0000000000,0.0002452658 +50.0000000000,2.0000000000,0.0004449008 +50.0000000000,3.0000000000,0.0004677162 +50.0000000000,4.0000000000,0.0006388319 +50.0000000000,5.0000000000,0.0008612822 +50.0000000000,6.0000000000,0.0008612822 +50.0000000000,7.0000000000,0.0009525439 +50.0000000000,8.0000000000,0.0009297285 +50.0000000000,9.0000000000,0.0010951403 +50.0000000000,10.0000000000,0.0011350673 +50.0000000000,11.0000000000,0.0009981748 +50.0000000000,12.0000000000,0.0008783938 +50.0000000000,13.0000000000,0.0007186858 +50.0000000000,14.0000000000,0.0004905316 +50.0000000000,15.0000000000,0.0003536391 +50.0000000000,16.0000000000,0.0001597080 +50.0000000000,17.0000000000,0.0000513347 +50.0000000000,18.0000000000,0.0000342231 +50.0000000000,19.0000000000,0.0000285193 +50.0000000000,20.0000000000,0.0000000000 +50.0000000000,21.0000000000,0.0000000000 +50.0000000000,22.0000000000,0.0000000000 +50.0000000000,23.0000000000,0.0000000000 +50.0000000000,24.0000000000,0.0000000000 +50.0000000000,25.0000000000,0.0000000000 +55.0000000000,0.0000000000,0.0000513347 +55.0000000000,1.0000000000,0.0001996350 +55.0000000000,2.0000000000,0.0003194159 +55.0000000000,3.0000000000,0.0004848277 +55.0000000000,4.0000000000,0.0006160164 +55.0000000000,5.0000000000,0.0007928360 +55.0000000000,6.0000000000,0.0008042437 +55.0000000000,7.0000000000,0.0010152863 +55.0000000000,8.0000000000,0.0009354324 +55.0000000000,9.0000000000,0.0009411362 +55.0000000000,10.0000000000,0.0009468401 +55.0000000000,11.0000000000,0.0009012092 +55.0000000000,12.0000000000,0.0007072781 +55.0000000000,13.0000000000,0.0005989049 +55.0000000000,14.0000000000,0.0004620123 +55.0000000000,15.0000000000,0.0004563085 +55.0000000000,16.0000000000,0.0002794889 +55.0000000000,17.0000000000,0.0000912617 +55.0000000000,18.0000000000,0.0000741501 +55.0000000000,19.0000000000,0.0000285193 +55.0000000000,20.0000000000,0.0000171116 +55.0000000000,21.0000000000,0.0000114077 +55.0000000000,22.0000000000,0.0000000000 +55.0000000000,23.0000000000,0.0000000000 +55.0000000000,24.0000000000,0.0000000000 +55.0000000000,25.0000000000,0.0000000000 +60.0000000000,0.0000000000,0.0000342231 +60.0000000000,1.0000000000,0.0002509697 +60.0000000000,2.0000000000,0.0003878622 +60.0000000000,3.0000000000,0.0005361624 +60.0000000000,4.0000000000,0.0006103126 +60.0000000000,5.0000000000,0.0008213552 +60.0000000000,6.0000000000,0.0008612822 +60.0000000000,7.0000000000,0.0008441707 +60.0000000000,8.0000000000,0.0010894365 +60.0000000000,9.0000000000,0.0009981748 +60.0000000000,10.0000000000,0.0011635866 +60.0000000000,11.0000000000,0.0010209902 +60.0000000000,12.0000000000,0.0009126169 +60.0000000000,13.0000000000,0.0008156514 +60.0000000000,14.0000000000,0.0004905316 +60.0000000000,15.0000000000,0.0004563085 +60.0000000000,16.0000000000,0.0002737851 +60.0000000000,17.0000000000,0.0001368925 +60.0000000000,18.0000000000,0.0000570386 +60.0000000000,19.0000000000,0.0000228154 +60.0000000000,20.0000000000,0.0000000000 +60.0000000000,21.0000000000,0.0000000000 +60.0000000000,22.0000000000,0.0000000000 +60.0000000000,23.0000000000,0.0000000000 +60.0000000000,24.0000000000,0.0000000000 +60.0000000000,25.0000000000,0.0000000000 +65.0000000000,0.0000000000,0.0000228154 +65.0000000000,1.0000000000,0.0001996350 +65.0000000000,2.0000000000,0.0002966005 +65.0000000000,3.0000000000,0.0004848277 +65.0000000000,4.0000000000,0.0007072781 +65.0000000000,5.0000000000,0.0007700205 +65.0000000000,6.0000000000,0.0009924709 +65.0000000000,7.0000000000,0.0010266940 +65.0000000000,8.0000000000,0.0010951403 +65.0000000000,9.0000000000,0.0011122519 +65.0000000000,10.0000000000,0.0010723249 +65.0000000000,11.0000000000,0.0010438056 +65.0000000000,12.0000000000,0.0008099475 +65.0000000000,13.0000000000,0.0007871321 +65.0000000000,14.0000000000,0.0004848277 +65.0000000000,15.0000000000,0.0004449008 +65.0000000000,16.0000000000,0.0001825234 +65.0000000000,17.0000000000,0.0001540041 +65.0000000000,18.0000000000,0.0000228154 +65.0000000000,19.0000000000,0.0000000000 +65.0000000000,20.0000000000,0.0000000000 +65.0000000000,21.0000000000,0.0000000000 +65.0000000000,22.0000000000,0.0000000000 +65.0000000000,23.0000000000,0.0000000000 +65.0000000000,24.0000000000,0.0000000000 +65.0000000000,25.0000000000,0.0000000000 +70.0000000000,0.0000000000,0.0000342231 +70.0000000000,1.0000000000,0.0001254848 +70.0000000000,2.0000000000,0.0002966005 +70.0000000000,3.0000000000,0.0005703856 +70.0000000000,4.0000000000,0.0005817933 +70.0000000000,5.0000000000,0.0008156514 +70.0000000000,6.0000000000,0.0009012092 +70.0000000000,7.0000000000,0.0012263290 +70.0000000000,8.0000000000,0.0010894365 +70.0000000000,9.0000000000,0.0012092174 +70.0000000000,10.0000000000,0.0010038786 +70.0000000000,11.0000000000,0.0012548483 +70.0000000000,12.0000000000,0.0009639516 +70.0000000000,13.0000000000,0.0005874971 +70.0000000000,14.0000000000,0.0004677162 +70.0000000000,15.0000000000,0.0003878622 +70.0000000000,16.0000000000,0.0001939311 +70.0000000000,17.0000000000,0.0001254848 +70.0000000000,18.0000000000,0.0000285193 +70.0000000000,19.0000000000,0.0000057039 +70.0000000000,20.0000000000,0.0000000000 +70.0000000000,21.0000000000,0.0000114077 +70.0000000000,22.0000000000,0.0000000000 +70.0000000000,23.0000000000,0.0000000000 +70.0000000000,24.0000000000,0.0000000000 +70.0000000000,25.0000000000,0.0000000000 +75.0000000000,0.0000000000,0.0000456308 +75.0000000000,1.0000000000,0.0002167465 +75.0000000000,2.0000000000,0.0003308236 +75.0000000000,3.0000000000,0.0004677162 +75.0000000000,4.0000000000,0.0006901666 +75.0000000000,5.0000000000,0.0006616473 +75.0000000000,6.0000000000,0.0009240246 +75.0000000000,7.0000000000,0.0010438056 +75.0000000000,8.0000000000,0.0011293634 +75.0000000000,9.0000000000,0.0012719598 +75.0000000000,10.0000000000,0.0012491444 +75.0000000000,11.0000000000,0.0010095825 +75.0000000000,12.0000000000,0.0008783938 +75.0000000000,13.0000000000,0.0006160164 +75.0000000000,14.0000000000,0.0004620123 +75.0000000000,15.0000000000,0.0003080082 +75.0000000000,16.0000000000,0.0001368925 +75.0000000000,17.0000000000,0.0001197810 +75.0000000000,18.0000000000,0.0000513347 +75.0000000000,19.0000000000,0.0000285193 +75.0000000000,20.0000000000,0.0000114077 +75.0000000000,21.0000000000,0.0000114077 +75.0000000000,22.0000000000,0.0000000000 +75.0000000000,23.0000000000,0.0000000000 +75.0000000000,24.0000000000,0.0000000000 +75.0000000000,25.0000000000,0.0000000000 +80.0000000000,0.0000000000,0.0000570386 +80.0000000000,1.0000000000,0.0001711157 +80.0000000000,2.0000000000,0.0003251198 +80.0000000000,3.0000000000,0.0003935661 +80.0000000000,4.0000000000,0.0005989049 +80.0000000000,5.0000000000,0.0008042437 +80.0000000000,6.0000000000,0.0009924709 +80.0000000000,7.0000000000,0.0012263290 +80.0000000000,8.0000000000,0.0011806982 +80.0000000000,9.0000000000,0.0012662560 +80.0000000000,10.0000000000,0.0010894365 +80.0000000000,11.0000000000,0.0011122519 +80.0000000000,12.0000000000,0.0008498745 +80.0000000000,13.0000000000,0.0006445357 +80.0000000000,14.0000000000,0.0005190509 +80.0000000000,15.0000000000,0.0003251198 +80.0000000000,16.0000000000,0.0001768195 +80.0000000000,17.0000000000,0.0001083733 +80.0000000000,18.0000000000,0.0000855578 +80.0000000000,19.0000000000,0.0000285193 +80.0000000000,20.0000000000,0.0000114077 +80.0000000000,21.0000000000,0.0000114077 +80.0000000000,22.0000000000,0.0000000000 +80.0000000000,23.0000000000,0.0000000000 +80.0000000000,24.0000000000,0.0000000000 +80.0000000000,25.0000000000,0.0000000000 +85.0000000000,0.0000000000,0.0000399270 +85.0000000000,1.0000000000,0.0002053388 +85.0000000000,2.0000000000,0.0002794889 +85.0000000000,3.0000000000,0.0005247547 +85.0000000000,4.0000000000,0.0005989049 +85.0000000000,5.0000000000,0.0007472051 +85.0000000000,6.0000000000,0.0010495095 +85.0000000000,7.0000000000,0.0012149213 +85.0000000000,8.0000000000,0.0012491444 +85.0000000000,9.0000000000,0.0013632215 +85.0000000000,10.0000000000,0.0012605521 +85.0000000000,11.0000000000,0.0010495095 +85.0000000000,12.0000000000,0.0007928360 +85.0000000000,13.0000000000,0.0006160164 +85.0000000000,14.0000000000,0.0004905316 +85.0000000000,15.0000000000,0.0003251198 +85.0000000000,16.0000000000,0.0002053388 +85.0000000000,17.0000000000,0.0001540041 +85.0000000000,18.0000000000,0.0000684463 +85.0000000000,19.0000000000,0.0000570386 +85.0000000000,20.0000000000,0.0000171116 +85.0000000000,21.0000000000,0.0000114077 +85.0000000000,22.0000000000,0.0000000000 +85.0000000000,23.0000000000,0.0000000000 +85.0000000000,24.0000000000,0.0000000000 +85.0000000000,25.0000000000,0.0000000000 +90.0000000000,0.0000000000,0.0000399270 +90.0000000000,1.0000000000,0.0001311887 +90.0000000000,2.0000000000,0.0002680812 +90.0000000000,3.0000000000,0.0004563085 +90.0000000000,4.0000000000,0.0005304586 +90.0000000000,5.0000000000,0.0008498745 +90.0000000000,6.0000000000,0.0008955054 +90.0000000000,7.0000000000,0.0010438056 +90.0000000000,8.0000000000,0.0011635866 +90.0000000000,9.0000000000,0.0011236596 +90.0000000000,10.0000000000,0.0010495095 +90.0000000000,11.0000000000,0.0008555784 +90.0000000000,12.0000000000,0.0007243897 +90.0000000000,13.0000000000,0.0006046087 +90.0000000000,14.0000000000,0.0003878622 +90.0000000000,15.0000000000,0.0002167465 +90.0000000000,16.0000000000,0.0002110427 +90.0000000000,17.0000000000,0.0001996350 +90.0000000000,18.0000000000,0.0000741501 +90.0000000000,19.0000000000,0.0000741501 +90.0000000000,20.0000000000,0.0000513347 +90.0000000000,21.0000000000,0.0000114077 +90.0000000000,22.0000000000,0.0000000000 +90.0000000000,23.0000000000,0.0000000000 +90.0000000000,24.0000000000,0.0000000000 +90.0000000000,25.0000000000,0.0000000000 +95.0000000000,0.0000000000,0.0000342231 +95.0000000000,1.0000000000,0.0002680812 +95.0000000000,2.0000000000,0.0003593429 +95.0000000000,3.0000000000,0.0005076432 +95.0000000000,4.0000000000,0.0005475702 +95.0000000000,5.0000000000,0.0005932010 +95.0000000000,6.0000000000,0.0009297285 +95.0000000000,7.0000000000,0.0009468401 +95.0000000000,8.0000000000,0.0010495095 +95.0000000000,9.0000000000,0.0010266940 +95.0000000000,10.0000000000,0.0009240246 +95.0000000000,11.0000000000,0.0009354324 +95.0000000000,12.0000000000,0.0008042437 +95.0000000000,13.0000000000,0.0006844627 +95.0000000000,14.0000000000,0.0005247547 +95.0000000000,15.0000000000,0.0004391969 +95.0000000000,16.0000000000,0.0002966005 +95.0000000000,17.0000000000,0.0001311887 +95.0000000000,18.0000000000,0.0000741501 +95.0000000000,19.0000000000,0.0000627424 +95.0000000000,20.0000000000,0.0000342231 +95.0000000000,21.0000000000,0.0000000000 +95.0000000000,22.0000000000,0.0000114077 +95.0000000000,23.0000000000,0.0000000000 +95.0000000000,24.0000000000,0.0000000000 +95.0000000000,25.0000000000,0.0000000000 +100.0000000000,0.0000000000,0.0000342231 +100.0000000000,1.0000000000,0.0001711157 +100.0000000000,2.0000000000,0.0003764545 +100.0000000000,3.0000000000,0.0004620123 +100.0000000000,4.0000000000,0.0005932010 +100.0000000000,5.0000000000,0.0007871321 +100.0000000000,6.0000000000,0.0008099475 +100.0000000000,7.0000000000,0.0009069131 +100.0000000000,8.0000000000,0.0008555784 +100.0000000000,9.0000000000,0.0009411362 +100.0000000000,10.0000000000,0.0008612822 +100.0000000000,11.0000000000,0.0008555784 +100.0000000000,12.0000000000,0.0007757244 +100.0000000000,13.0000000000,0.0007357974 +100.0000000000,14.0000000000,0.0004734200 +100.0000000000,15.0000000000,0.0004049738 +100.0000000000,16.0000000000,0.0002053388 +100.0000000000,17.0000000000,0.0000741501 +100.0000000000,18.0000000000,0.0000513347 +100.0000000000,19.0000000000,0.0000114077 +100.0000000000,20.0000000000,0.0000057039 +100.0000000000,21.0000000000,0.0000000000 +100.0000000000,22.0000000000,0.0000000000 +100.0000000000,23.0000000000,0.0000000000 +100.0000000000,24.0000000000,0.0000000000 +100.0000000000,25.0000000000,0.0000000000 +105.0000000000,0.0000000000,0.0000114077 +105.0000000000,1.0000000000,0.0002509697 +105.0000000000,2.0000000000,0.0003536391 +105.0000000000,3.0000000000,0.0004848277 +105.0000000000,4.0000000000,0.0006616473 +105.0000000000,5.0000000000,0.0006616473 +105.0000000000,6.0000000000,0.0007928360 +105.0000000000,7.0000000000,0.0007985398 +105.0000000000,8.0000000000,0.0008042437 +105.0000000000,9.0000000000,0.0009069131 +105.0000000000,10.0000000000,0.0007757244 +105.0000000000,11.0000000000,0.0007072781 +105.0000000000,12.0000000000,0.0006331280 +105.0000000000,13.0000000000,0.0003935661 +105.0000000000,14.0000000000,0.0004049738 +105.0000000000,15.0000000000,0.0002338581 +105.0000000000,16.0000000000,0.0001996350 +105.0000000000,17.0000000000,0.0000741501 +105.0000000000,18.0000000000,0.0000228154 +105.0000000000,19.0000000000,0.0000000000 +105.0000000000,20.0000000000,0.0000000000 +105.0000000000,21.0000000000,0.0000000000 +105.0000000000,22.0000000000,0.0000057039 +105.0000000000,23.0000000000,0.0000000000 +105.0000000000,24.0000000000,0.0000000000 +105.0000000000,25.0000000000,0.0000000000 +110.0000000000,0.0000000000,0.0000342231 +110.0000000000,1.0000000000,0.0002110427 +110.0000000000,2.0000000000,0.0003479352 +110.0000000000,3.0000000000,0.0004391969 +110.0000000000,4.0000000000,0.0005475702 +110.0000000000,5.0000000000,0.0007072781 +110.0000000000,6.0000000000,0.0007015743 +110.0000000000,7.0000000000,0.0006901666 +110.0000000000,8.0000000000,0.0007928360 +110.0000000000,9.0000000000,0.0008726899 +110.0000000000,10.0000000000,0.0006958704 +110.0000000000,11.0000000000,0.0005703856 +110.0000000000,12.0000000000,0.0005019393 +110.0000000000,13.0000000000,0.0002966005 +110.0000000000,14.0000000000,0.0002908966 +110.0000000000,15.0000000000,0.0001654118 +110.0000000000,16.0000000000,0.0001483003 +110.0000000000,17.0000000000,0.0000399270 +110.0000000000,18.0000000000,0.0000456308 +110.0000000000,19.0000000000,0.0000057039 +110.0000000000,20.0000000000,0.0000057039 +110.0000000000,21.0000000000,0.0000000000 +110.0000000000,22.0000000000,0.0000000000 +110.0000000000,23.0000000000,0.0000000000 +110.0000000000,24.0000000000,0.0000000000 +110.0000000000,25.0000000000,0.0000000000 +115.0000000000,0.0000000000,0.0000285193 +115.0000000000,1.0000000000,0.0001425964 +115.0000000000,2.0000000000,0.0003422313 +115.0000000000,3.0000000000,0.0004506046 +115.0000000000,4.0000000000,0.0004905316 +115.0000000000,5.0000000000,0.0005361624 +115.0000000000,6.0000000000,0.0006502396 +115.0000000000,7.0000000000,0.0007928360 +115.0000000000,8.0000000000,0.0007643167 +115.0000000000,9.0000000000,0.0007757244 +115.0000000000,10.0000000000,0.0006445357 +115.0000000000,11.0000000000,0.0006673511 +115.0000000000,12.0000000000,0.0005646817 +115.0000000000,13.0000000000,0.0004106776 +115.0000000000,14.0000000000,0.0003479352 +115.0000000000,15.0000000000,0.0002509697 +115.0000000000,16.0000000000,0.0001825234 +115.0000000000,17.0000000000,0.0000513347 +115.0000000000,18.0000000000,0.0000228154 +115.0000000000,19.0000000000,0.0000171116 +115.0000000000,20.0000000000,0.0000000000 +115.0000000000,21.0000000000,0.0000000000 +115.0000000000,22.0000000000,0.0000000000 +115.0000000000,23.0000000000,0.0000000000 +115.0000000000,24.0000000000,0.0000000000 +115.0000000000,25.0000000000,0.0000000000 +120.0000000000,0.0000000000,0.0000171116 +120.0000000000,1.0000000000,0.0001311887 +120.0000000000,2.0000000000,0.0003422313 +120.0000000000,3.0000000000,0.0004049738 +120.0000000000,4.0000000000,0.0005589779 +120.0000000000,5.0000000000,0.0006046087 +120.0000000000,6.0000000000,0.0008099475 +120.0000000000,7.0000000000,0.0008783938 +120.0000000000,8.0000000000,0.0007985398 +120.0000000000,9.0000000000,0.0008955054 +120.0000000000,10.0000000000,0.0007700205 +120.0000000000,11.0000000000,0.0006388319 +120.0000000000,12.0000000000,0.0005475702 +120.0000000000,13.0000000000,0.0003593429 +120.0000000000,14.0000000000,0.0002851928 +120.0000000000,15.0000000000,0.0002908966 +120.0000000000,16.0000000000,0.0001939311 +120.0000000000,17.0000000000,0.0000741501 +120.0000000000,18.0000000000,0.0000399270 +120.0000000000,19.0000000000,0.0000057039 +120.0000000000,20.0000000000,0.0000057039 +120.0000000000,21.0000000000,0.0000000000 +120.0000000000,22.0000000000,0.0000000000 +120.0000000000,23.0000000000,0.0000000000 +120.0000000000,24.0000000000,0.0000000000 +120.0000000000,25.0000000000,0.0000000000 +125.0000000000,0.0000000000,0.0000399270 +125.0000000000,1.0000000000,0.0001483003 +125.0000000000,2.0000000000,0.0003023044 +125.0000000000,3.0000000000,0.0004563085 +125.0000000000,4.0000000000,0.0007129820 +125.0000000000,5.0000000000,0.0006673511 +125.0000000000,6.0000000000,0.0007586128 +125.0000000000,7.0000000000,0.0006901666 +125.0000000000,8.0000000000,0.0007757244 +125.0000000000,9.0000000000,0.0007529090 +125.0000000000,10.0000000000,0.0008213552 +125.0000000000,11.0000000000,0.0007015743 +125.0000000000,12.0000000000,0.0004734200 +125.0000000000,13.0000000000,0.0004106776 +125.0000000000,14.0000000000,0.0002680812 +125.0000000000,15.0000000000,0.0002110427 +125.0000000000,16.0000000000,0.0001026694 +125.0000000000,17.0000000000,0.0000627424 +125.0000000000,18.0000000000,0.0000342231 +125.0000000000,19.0000000000,0.0000228154 +125.0000000000,20.0000000000,0.0000057039 +125.0000000000,21.0000000000,0.0000057039 +125.0000000000,22.0000000000,0.0000000000 +125.0000000000,23.0000000000,0.0000000000 +125.0000000000,24.0000000000,0.0000000000 +125.0000000000,25.0000000000,0.0000000000 +130.0000000000,0.0000000000,0.0000285193 +130.0000000000,1.0000000000,0.0001768195 +130.0000000000,2.0000000000,0.0002509697 +130.0000000000,3.0000000000,0.0005361624 +130.0000000000,4.0000000000,0.0006388319 +130.0000000000,5.0000000000,0.0005589779 +130.0000000000,6.0000000000,0.0007814282 +130.0000000000,7.0000000000,0.0007529090 +130.0000000000,8.0000000000,0.0008726899 +130.0000000000,9.0000000000,0.0007357974 +130.0000000000,10.0000000000,0.0006274241 +130.0000000000,11.0000000000,0.0006502396 +130.0000000000,12.0000000000,0.0005019393 +130.0000000000,13.0000000000,0.0004277892 +130.0000000000,14.0000000000,0.0002966005 +130.0000000000,15.0000000000,0.0002224504 +130.0000000000,16.0000000000,0.0001483003 +130.0000000000,17.0000000000,0.0000570386 +130.0000000000,18.0000000000,0.0000114077 +130.0000000000,19.0000000000,0.0000114077 +130.0000000000,20.0000000000,0.0000114077 +130.0000000000,21.0000000000,0.0000000000 +130.0000000000,22.0000000000,0.0000000000 +130.0000000000,23.0000000000,0.0000000000 +130.0000000000,24.0000000000,0.0000000000 +130.0000000000,25.0000000000,0.0000000000 +135.0000000000,0.0000000000,0.0000684463 +135.0000000000,1.0000000000,0.0001768195 +135.0000000000,2.0000000000,0.0002908966 +135.0000000000,3.0000000000,0.0005019393 +135.0000000000,4.0000000000,0.0005760894 +135.0000000000,5.0000000000,0.0005874971 +135.0000000000,6.0000000000,0.0006673511 +135.0000000000,7.0000000000,0.0007814282 +135.0000000000,8.0000000000,0.0008042437 +135.0000000000,9.0000000000,0.0007357974 +135.0000000000,10.0000000000,0.0005989049 +135.0000000000,11.0000000000,0.0005589779 +135.0000000000,12.0000000000,0.0003479352 +135.0000000000,13.0000000000,0.0004049738 +135.0000000000,14.0000000000,0.0002623774 +135.0000000000,15.0000000000,0.0001425964 +135.0000000000,16.0000000000,0.0001368925 +135.0000000000,17.0000000000,0.0000513347 +135.0000000000,18.0000000000,0.0000285193 +135.0000000000,19.0000000000,0.0000114077 +135.0000000000,20.0000000000,0.0000000000 +135.0000000000,21.0000000000,0.0000000000 +135.0000000000,22.0000000000,0.0000000000 +135.0000000000,23.0000000000,0.0000000000 +135.0000000000,24.0000000000,0.0000000000 +135.0000000000,25.0000000000,0.0000000000 +140.0000000000,0.0000000000,0.0000342231 +140.0000000000,1.0000000000,0.0001425964 +140.0000000000,2.0000000000,0.0003137121 +140.0000000000,3.0000000000,0.0005589779 +140.0000000000,4.0000000000,0.0005760894 +140.0000000000,5.0000000000,0.0006331280 +140.0000000000,6.0000000000,0.0007129820 +140.0000000000,7.0000000000,0.0006388319 +140.0000000000,8.0000000000,0.0006616473 +140.0000000000,9.0000000000,0.0006559434 +140.0000000000,10.0000000000,0.0005932010 +140.0000000000,11.0000000000,0.0005361624 +140.0000000000,12.0000000000,0.0004848277 +140.0000000000,13.0000000000,0.0003593429 +140.0000000000,14.0000000000,0.0002623774 +140.0000000000,15.0000000000,0.0001083733 +140.0000000000,16.0000000000,0.0000855578 +140.0000000000,17.0000000000,0.0000513347 +140.0000000000,18.0000000000,0.0000342231 +140.0000000000,19.0000000000,0.0000114077 +140.0000000000,20.0000000000,0.0000114077 +140.0000000000,21.0000000000,0.0000000000 +140.0000000000,22.0000000000,0.0000000000 +140.0000000000,23.0000000000,0.0000000000 +140.0000000000,24.0000000000,0.0000000000 +140.0000000000,25.0000000000,0.0000000000 +145.0000000000,0.0000000000,0.0000342231 +145.0000000000,1.0000000000,0.0002224504 +145.0000000000,2.0000000000,0.0002566735 +145.0000000000,3.0000000000,0.0005190509 +145.0000000000,4.0000000000,0.0005589779 +145.0000000000,5.0000000000,0.0006046087 +145.0000000000,6.0000000000,0.0007586128 +145.0000000000,7.0000000000,0.0006274241 +145.0000000000,8.0000000000,0.0006901666 +145.0000000000,9.0000000000,0.0006217203 +145.0000000000,10.0000000000,0.0005932010 +145.0000000000,11.0000000000,0.0005475702 +145.0000000000,12.0000000000,0.0004848277 +145.0000000000,13.0000000000,0.0004049738 +145.0000000000,14.0000000000,0.0002281542 +145.0000000000,15.0000000000,0.0001026694 +145.0000000000,16.0000000000,0.0000912617 +145.0000000000,17.0000000000,0.0000570386 +145.0000000000,18.0000000000,0.0000285193 +145.0000000000,19.0000000000,0.0000228154 +145.0000000000,20.0000000000,0.0000114077 +145.0000000000,21.0000000000,0.0000057039 +145.0000000000,22.0000000000,0.0000000000 +145.0000000000,23.0000000000,0.0000000000 +145.0000000000,24.0000000000,0.0000000000 +145.0000000000,25.0000000000,0.0000000000 +150.0000000000,0.0000000000,0.0000285193 +150.0000000000,1.0000000000,0.0001768195 +150.0000000000,2.0000000000,0.0003536391 +150.0000000000,3.0000000000,0.0005989049 +150.0000000000,4.0000000000,0.0005361624 +150.0000000000,5.0000000000,0.0006160164 +150.0000000000,6.0000000000,0.0006844627 +150.0000000000,7.0000000000,0.0006730550 +150.0000000000,8.0000000000,0.0006388319 +150.0000000000,9.0000000000,0.0006787588 +150.0000000000,10.0000000000,0.0007072781 +150.0000000000,11.0000000000,0.0006844627 +150.0000000000,12.0000000000,0.0004734200 +150.0000000000,13.0000000000,0.0003080082 +150.0000000000,14.0000000000,0.0002623774 +150.0000000000,15.0000000000,0.0001540041 +150.0000000000,16.0000000000,0.0001654118 +150.0000000000,17.0000000000,0.0000798540 +150.0000000000,18.0000000000,0.0000513347 +150.0000000000,19.0000000000,0.0000171116 +150.0000000000,20.0000000000,0.0000228154 +150.0000000000,21.0000000000,0.0000057039 +150.0000000000,22.0000000000,0.0000057039 +150.0000000000,23.0000000000,0.0000000000 +150.0000000000,24.0000000000,0.0000000000 +150.0000000000,25.0000000000,0.0000000000 +155.0000000000,0.0000000000,0.0000342231 +155.0000000000,1.0000000000,0.0001996350 +155.0000000000,2.0000000000,0.0002966005 +155.0000000000,3.0000000000,0.0005133470 +155.0000000000,4.0000000000,0.0005646817 +155.0000000000,5.0000000000,0.0007186858 +155.0000000000,6.0000000000,0.0006445357 +155.0000000000,7.0000000000,0.0007357974 +155.0000000000,8.0000000000,0.0005932010 +155.0000000000,9.0000000000,0.0007928360 +155.0000000000,10.0000000000,0.0007243897 +155.0000000000,11.0000000000,0.0005760894 +155.0000000000,12.0000000000,0.0004506046 +155.0000000000,13.0000000000,0.0003650468 +155.0000000000,14.0000000000,0.0003251198 +155.0000000000,15.0000000000,0.0002680812 +155.0000000000,16.0000000000,0.0001140771 +155.0000000000,17.0000000000,0.0000855578 +155.0000000000,18.0000000000,0.0000513347 +155.0000000000,19.0000000000,0.0000171116 +155.0000000000,20.0000000000,0.0000342231 +155.0000000000,21.0000000000,0.0000114077 +155.0000000000,22.0000000000,0.0000114077 +155.0000000000,23.0000000000,0.0000000000 +155.0000000000,24.0000000000,0.0000057039 +155.0000000000,25.0000000000,0.0000000000 +160.0000000000,0.0000000000,0.0000399270 +160.0000000000,1.0000000000,0.0001197810 +160.0000000000,2.0000000000,0.0003650468 +160.0000000000,3.0000000000,0.0003935661 +160.0000000000,4.0000000000,0.0005475702 +160.0000000000,5.0000000000,0.0007415013 +160.0000000000,6.0000000000,0.0006217203 +160.0000000000,7.0000000000,0.0007015743 +160.0000000000,8.0000000000,0.0007529090 +160.0000000000,9.0000000000,0.0007357974 +160.0000000000,10.0000000000,0.0007186858 +160.0000000000,11.0000000000,0.0006787588 +160.0000000000,12.0000000000,0.0005646817 +160.0000000000,13.0000000000,0.0004905316 +160.0000000000,14.0000000000,0.0003764545 +160.0000000000,15.0000000000,0.0003080082 +160.0000000000,16.0000000000,0.0001825234 +160.0000000000,17.0000000000,0.0001311887 +160.0000000000,18.0000000000,0.0000969655 +160.0000000000,19.0000000000,0.0000684463 +160.0000000000,20.0000000000,0.0000171116 +160.0000000000,21.0000000000,0.0000171116 +160.0000000000,22.0000000000,0.0000000000 +160.0000000000,23.0000000000,0.0000000000 +160.0000000000,24.0000000000,0.0000057039 +160.0000000000,25.0000000000,0.0000000000 +165.0000000000,0.0000000000,0.0000399270 +165.0000000000,1.0000000000,0.0001825234 +165.0000000000,2.0000000000,0.0003878622 +165.0000000000,3.0000000000,0.0005361624 +165.0000000000,4.0000000000,0.0005190509 +165.0000000000,5.0000000000,0.0006103126 +165.0000000000,6.0000000000,0.0007243897 +165.0000000000,7.0000000000,0.0006331280 +165.0000000000,8.0000000000,0.0007586128 +165.0000000000,9.0000000000,0.0007472051 +165.0000000000,10.0000000000,0.0007300935 +165.0000000000,11.0000000000,0.0007871321 +165.0000000000,12.0000000000,0.0006502396 +165.0000000000,13.0000000000,0.0005703856 +165.0000000000,14.0000000000,0.0004962355 +165.0000000000,15.0000000000,0.0003194159 +165.0000000000,16.0000000000,0.0003080082 +165.0000000000,17.0000000000,0.0002053388 +165.0000000000,18.0000000000,0.0001425964 +165.0000000000,19.0000000000,0.0001197810 +165.0000000000,20.0000000000,0.0000513347 +165.0000000000,21.0000000000,0.0000342231 +165.0000000000,22.0000000000,0.0000057039 +165.0000000000,23.0000000000,0.0000000000 +165.0000000000,24.0000000000,0.0000000000 +165.0000000000,25.0000000000,0.0000000000 +170.0000000000,0.0000000000,0.0000513347 +170.0000000000,1.0000000000,0.0001311887 +170.0000000000,2.0000000000,0.0002680812 +170.0000000000,3.0000000000,0.0004734200 +170.0000000000,4.0000000000,0.0006160164 +170.0000000000,5.0000000000,0.0005989049 +170.0000000000,6.0000000000,0.0007586128 +170.0000000000,7.0000000000,0.0009012092 +170.0000000000,8.0000000000,0.0008384668 +170.0000000000,9.0000000000,0.0009411362 +170.0000000000,10.0000000000,0.0008555784 +170.0000000000,11.0000000000,0.0009126169 +170.0000000000,12.0000000000,0.0007871321 +170.0000000000,13.0000000000,0.0006559434 +170.0000000000,14.0000000000,0.0005989049 +170.0000000000,15.0000000000,0.0004962355 +170.0000000000,16.0000000000,0.0004391969 +170.0000000000,17.0000000000,0.0003536391 +170.0000000000,18.0000000000,0.0002053388 +170.0000000000,19.0000000000,0.0001425964 +170.0000000000,20.0000000000,0.0001026694 +170.0000000000,21.0000000000,0.0000570386 +170.0000000000,22.0000000000,0.0000171116 +170.0000000000,23.0000000000,0.0000057039 +170.0000000000,24.0000000000,0.0000057039 +170.0000000000,25.0000000000,0.0000057039 +175.0000000000,0.0000000000,0.0000285193 +175.0000000000,1.0000000000,0.0001711157 +175.0000000000,2.0000000000,0.0003023044 +175.0000000000,3.0000000000,0.0004220853 +175.0000000000,4.0000000000,0.0005475702 +175.0000000000,5.0000000000,0.0005532740 +175.0000000000,6.0000000000,0.0007814282 +175.0000000000,7.0000000000,0.0008783938 +175.0000000000,8.0000000000,0.0008840977 +175.0000000000,9.0000000000,0.0008840977 +175.0000000000,10.0000000000,0.0010495095 +175.0000000000,11.0000000000,0.0008840977 +175.0000000000,12.0000000000,0.0007472051 +175.0000000000,13.0000000000,0.0006046087 +175.0000000000,14.0000000000,0.0007072781 +175.0000000000,15.0000000000,0.0005475702 +175.0000000000,16.0000000000,0.0003878622 +175.0000000000,17.0000000000,0.0003650468 +175.0000000000,18.0000000000,0.0003194159 +175.0000000000,19.0000000000,0.0002053388 +175.0000000000,20.0000000000,0.0000912617 +175.0000000000,21.0000000000,0.0000912617 +175.0000000000,22.0000000000,0.0000228154 +175.0000000000,23.0000000000,0.0000114077 +175.0000000000,24.0000000000,0.0000057039 +175.0000000000,25.0000000000,0.0000171116 +180.0000000000,0.0000000000,0.0000285193 +180.0000000000,1.0000000000,0.0001597080 +180.0000000000,2.0000000000,0.0002737851 +180.0000000000,3.0000000000,0.0005247547 +180.0000000000,4.0000000000,0.0005532740 +180.0000000000,5.0000000000,0.0008612822 +180.0000000000,6.0000000000,0.0007129820 +180.0000000000,7.0000000000,0.0008156514 +180.0000000000,8.0000000000,0.0009354324 +180.0000000000,9.0000000000,0.0009411362 +180.0000000000,10.0000000000,0.0008898015 +180.0000000000,11.0000000000,0.0008213552 +180.0000000000,12.0000000000,0.0009696555 +180.0000000000,13.0000000000,0.0007700205 +180.0000000000,14.0000000000,0.0006787588 +180.0000000000,15.0000000000,0.0006160164 +180.0000000000,16.0000000000,0.0004620123 +180.0000000000,17.0000000000,0.0004220853 +180.0000000000,18.0000000000,0.0001654118 +180.0000000000,19.0000000000,0.0001939311 +180.0000000000,20.0000000000,0.0001768195 +180.0000000000,21.0000000000,0.0001311887 +180.0000000000,22.0000000000,0.0000456308 +180.0000000000,23.0000000000,0.0000171116 +180.0000000000,24.0000000000,0.0000114077 +180.0000000000,25.0000000000,0.0000228154 +185.0000000000,0.0000000000,0.0000228154 +185.0000000000,1.0000000000,0.0001540041 +185.0000000000,2.0000000000,0.0002737851 +185.0000000000,3.0000000000,0.0004449008 +185.0000000000,4.0000000000,0.0004449008 +185.0000000000,5.0000000000,0.0009126169 +185.0000000000,6.0000000000,0.0008213552 +185.0000000000,7.0000000000,0.0011008442 +185.0000000000,8.0000000000,0.0010666210 +185.0000000000,9.0000000000,0.0010951403 +185.0000000000,10.0000000000,0.0009696555 +185.0000000000,11.0000000000,0.0009411362 +185.0000000000,12.0000000000,0.0009582478 +185.0000000000,13.0000000000,0.0008156514 +185.0000000000,14.0000000000,0.0007814282 +185.0000000000,15.0000000000,0.0007814282 +185.0000000000,16.0000000000,0.0006445357 +185.0000000000,17.0000000000,0.0004677162 +185.0000000000,18.0000000000,0.0003707506 +185.0000000000,19.0000000000,0.0003650468 +185.0000000000,20.0000000000,0.0002566735 +185.0000000000,21.0000000000,0.0001825234 +185.0000000000,22.0000000000,0.0000741501 +185.0000000000,23.0000000000,0.0000342231 +185.0000000000,24.0000000000,0.0000171116 +185.0000000000,25.0000000000,0.0000171116 +190.0000000000,0.0000000000,0.0000456308 +190.0000000000,1.0000000000,0.0001425964 +190.0000000000,2.0000000000,0.0003080082 +190.0000000000,3.0000000000,0.0003764545 +190.0000000000,4.0000000000,0.0004962355 +190.0000000000,5.0000000000,0.0005932010 +190.0000000000,6.0000000000,0.0009240246 +190.0000000000,7.0000000000,0.0010323979 +190.0000000000,8.0000000000,0.0010837326 +190.0000000000,9.0000000000,0.0010438056 +190.0000000000,10.0000000000,0.0011293634 +190.0000000000,11.0000000000,0.0010552133 +190.0000000000,12.0000000000,0.0009126169 +190.0000000000,13.0000000000,0.0009981748 +190.0000000000,14.0000000000,0.0009297285 +190.0000000000,15.0000000000,0.0007472051 +190.0000000000,16.0000000000,0.0007015743 +190.0000000000,17.0000000000,0.0007072781 +190.0000000000,18.0000000000,0.0005817933 +190.0000000000,19.0000000000,0.0004791239 +190.0000000000,20.0000000000,0.0002509697 +190.0000000000,21.0000000000,0.0001711157 +190.0000000000,22.0000000000,0.0001083733 +190.0000000000,23.0000000000,0.0000456308 +190.0000000000,24.0000000000,0.0000285193 +190.0000000000,25.0000000000,0.0000171116 +195.0000000000,0.0000000000,0.0000342231 +195.0000000000,1.0000000000,0.0001825234 +195.0000000000,2.0000000000,0.0003023044 +195.0000000000,3.0000000000,0.0004563085 +195.0000000000,4.0000000000,0.0005989049 +195.0000000000,5.0000000000,0.0006787588 +195.0000000000,6.0000000000,0.0008840977 +195.0000000000,7.0000000000,0.0008898015 +195.0000000000,8.0000000000,0.0010837326 +195.0000000000,9.0000000000,0.0010038786 +195.0000000000,10.0000000000,0.0011407712 +195.0000000000,11.0000000000,0.0012548483 +195.0000000000,12.0000000000,0.0010323979 +195.0000000000,13.0000000000,0.0009696555 +195.0000000000,14.0000000000,0.0009468401 +195.0000000000,15.0000000000,0.0008099475 +195.0000000000,16.0000000000,0.0008498745 +195.0000000000,17.0000000000,0.0006274241 +195.0000000000,18.0000000000,0.0004506046 +195.0000000000,19.0000000000,0.0004791239 +195.0000000000,20.0000000000,0.0002851928 +195.0000000000,21.0000000000,0.0002110427 +195.0000000000,22.0000000000,0.0001197810 +195.0000000000,23.0000000000,0.0000513347 +195.0000000000,24.0000000000,0.0000627424 +195.0000000000,25.0000000000,0.0000342231 +200.0000000000,0.0000000000,0.0000285193 +200.0000000000,1.0000000000,0.0001368925 +200.0000000000,2.0000000000,0.0003764545 +200.0000000000,3.0000000000,0.0005247547 +200.0000000000,4.0000000000,0.0006958704 +200.0000000000,5.0000000000,0.0007015743 +200.0000000000,6.0000000000,0.0008726899 +200.0000000000,7.0000000000,0.0011635866 +200.0000000000,8.0000000000,0.0011464750 +200.0000000000,9.0000000000,0.0012263290 +200.0000000000,10.0000000000,0.0011521789 +200.0000000000,11.0000000000,0.0012320329 +200.0000000000,12.0000000000,0.0012263290 +200.0000000000,13.0000000000,0.0011407712 +200.0000000000,14.0000000000,0.0012377367 +200.0000000000,15.0000000000,0.0010894365 +200.0000000000,16.0000000000,0.0010152863 +200.0000000000,17.0000000000,0.0008669861 +200.0000000000,18.0000000000,0.0007814282 +200.0000000000,19.0000000000,0.0006103126 +200.0000000000,20.0000000000,0.0004277892 +200.0000000000,21.0000000000,0.0003365275 +200.0000000000,22.0000000000,0.0002452658 +200.0000000000,23.0000000000,0.0001425964 +200.0000000000,24.0000000000,0.0000342231 +200.0000000000,25.0000000000,0.0000285193 +205.0000000000,0.0000000000,0.0000114077 +205.0000000000,1.0000000000,0.0001540041 +205.0000000000,2.0000000000,0.0003878622 +205.0000000000,3.0000000000,0.0005361624 +205.0000000000,4.0000000000,0.0006673511 +205.0000000000,5.0000000000,0.0008270591 +205.0000000000,6.0000000000,0.0009525439 +205.0000000000,7.0000000000,0.0011578827 +205.0000000000,8.0000000000,0.0013746292 +205.0000000000,9.0000000000,0.0015457449 +205.0000000000,10.0000000000,0.0013632215 +205.0000000000,11.0000000000,0.0014658909 +205.0000000000,12.0000000000,0.0012890714 +205.0000000000,13.0000000000,0.0013461100 +205.0000000000,14.0000000000,0.0011806982 +205.0000000000,15.0000000000,0.0010780287 +205.0000000000,16.0000000000,0.0009126169 +205.0000000000,17.0000000000,0.0010951403 +205.0000000000,18.0000000000,0.0008726899 +205.0000000000,19.0000000000,0.0006388319 +205.0000000000,20.0000000000,0.0005361624 +205.0000000000,21.0000000000,0.0003308236 +205.0000000000,22.0000000000,0.0001939311 +205.0000000000,23.0000000000,0.0001254848 +205.0000000000,24.0000000000,0.0000741501 +205.0000000000,25.0000000000,0.0000456308 +210.0000000000,0.0000000000,0.0000456308 +210.0000000000,1.0000000000,0.0001483003 +210.0000000000,2.0000000000,0.0003308236 +210.0000000000,3.0000000000,0.0006331280 +210.0000000000,4.0000000000,0.0006844627 +210.0000000000,5.0000000000,0.0009696555 +210.0000000000,6.0000000000,0.0013118868 +210.0000000000,7.0000000000,0.0012320329 +210.0000000000,8.0000000000,0.0013061830 +210.0000000000,9.0000000000,0.0015343372 +210.0000000000,10.0000000000,0.0017111567 +210.0000000000,11.0000000000,0.0016255989 +210.0000000000,12.0000000000,0.0016141912 +210.0000000000,13.0000000000,0.0015685603 +210.0000000000,14.0000000000,0.0014658909 +210.0000000000,15.0000000000,0.0013689254 +210.0000000000,16.0000000000,0.0012662560 +210.0000000000,17.0000000000,0.0012662560 +210.0000000000,18.0000000000,0.0009411362 +210.0000000000,19.0000000000,0.0007300935 +210.0000000000,20.0000000000,0.0006103126 +210.0000000000,21.0000000000,0.0003479352 +210.0000000000,22.0000000000,0.0002851928 +210.0000000000,23.0000000000,0.0001597080 +210.0000000000,24.0000000000,0.0001026694 +210.0000000000,25.0000000000,0.0000399270 +215.0000000000,0.0000000000,0.0000456308 +215.0000000000,1.0000000000,0.0001825234 +215.0000000000,2.0000000000,0.0004106776 +215.0000000000,3.0000000000,0.0004734200 +215.0000000000,4.0000000000,0.0008099475 +215.0000000000,5.0000000000,0.0010266940 +215.0000000000,6.0000000000,0.0011692904 +215.0000000000,7.0000000000,0.0013575177 +215.0000000000,8.0000000000,0.0016598220 +215.0000000000,9.0000000000,0.0017510837 +215.0000000000,10.0000000000,0.0021389459 +215.0000000000,11.0000000000,0.0019221994 +215.0000000000,12.0000000000,0.0020077572 +215.0000000000,13.0000000000,0.0020305727 +215.0000000000,14.0000000000,0.0016826375 +215.0000000000,15.0000000000,0.0015856719 +215.0000000000,16.0000000000,0.0015457449 +215.0000000000,17.0000000000,0.0012890714 +215.0000000000,18.0000000000,0.0010552133 +215.0000000000,19.0000000000,0.0008156514 +215.0000000000,20.0000000000,0.0004962355 +215.0000000000,21.0000000000,0.0003992699 +215.0000000000,22.0000000000,0.0002908966 +215.0000000000,23.0000000000,0.0002623774 +215.0000000000,24.0000000000,0.0001140771 +215.0000000000,25.0000000000,0.0000570386 +220.0000000000,0.0000000000,0.0000456308 +220.0000000000,1.0000000000,0.0001996350 +220.0000000000,2.0000000000,0.0003878622 +220.0000000000,3.0000000000,0.0006046087 +220.0000000000,4.0000000000,0.0010095825 +220.0000000000,5.0000000000,0.0010266940 +220.0000000000,6.0000000000,0.0012719598 +220.0000000000,7.0000000000,0.0013289984 +220.0000000000,8.0000000000,0.0017111567 +220.0000000000,9.0000000000,0.0020647958 +220.0000000000,10.0000000000,0.0021902806 +220.0000000000,11.0000000000,0.0024412503 +220.0000000000,12.0000000000,0.0024013233 +220.0000000000,13.0000000000,0.0022073922 +220.0000000000,14.0000000000,0.0019050878 +220.0000000000,15.0000000000,0.0020362765 +220.0000000000,16.0000000000,0.0017453799 +220.0000000000,17.0000000000,0.0013860370 +220.0000000000,18.0000000000,0.0011008442 +220.0000000000,19.0000000000,0.0008384668 +220.0000000000,20.0000000000,0.0007472051 +220.0000000000,21.0000000000,0.0005989049 +220.0000000000,22.0000000000,0.0003878622 +220.0000000000,23.0000000000,0.0002737851 +220.0000000000,24.0000000000,0.0001311887 +220.0000000000,25.0000000000,0.0000627424 +225.0000000000,0.0000000000,0.0000399270 +225.0000000000,1.0000000000,0.0001882272 +225.0000000000,2.0000000000,0.0003365275 +225.0000000000,3.0000000000,0.0006388319 +225.0000000000,4.0000000000,0.0008327629 +225.0000000000,5.0000000000,0.0009069131 +225.0000000000,6.0000000000,0.0015115218 +225.0000000000,7.0000000000,0.0015343372 +225.0000000000,8.0000000000,0.0019450148 +225.0000000000,9.0000000000,0.0024013233 +225.0000000000,10.0000000000,0.0024184349 +225.0000000000,11.0000000000,0.0025781428 +225.0000000000,12.0000000000,0.0027264431 +225.0000000000,13.0000000000,0.0024412503 +225.0000000000,14.0000000000,0.0023613963 +225.0000000000,15.0000000000,0.0019906457 +225.0000000000,16.0000000000,0.0019621264 +225.0000000000,17.0000000000,0.0014316678 +225.0000000000,18.0000000000,0.0011008442 +225.0000000000,19.0000000000,0.0009069131 +225.0000000000,20.0000000000,0.0007415013 +225.0000000000,21.0000000000,0.0005817933 +225.0000000000,22.0000000000,0.0003821583 +225.0000000000,23.0000000000,0.0002110427 +225.0000000000,24.0000000000,0.0001311887 +225.0000000000,25.0000000000,0.0000627424 +230.0000000000,0.0000000000,0.0000399270 +230.0000000000,1.0000000000,0.0002395619 +230.0000000000,2.0000000000,0.0003650468 +230.0000000000,3.0000000000,0.0006046087 +230.0000000000,4.0000000000,0.0009354324 +230.0000000000,5.0000000000,0.0011578827 +230.0000000000,6.0000000000,0.0012719598 +230.0000000000,7.0000000000,0.0017453799 +230.0000000000,8.0000000000,0.0018252339 +230.0000000000,9.0000000000,0.0021731691 +230.0000000000,10.0000000000,0.0026009582 +230.0000000000,11.0000000000,0.0025268081 +230.0000000000,12.0000000000,0.0023842117 +230.0000000000,13.0000000000,0.0022986539 +230.0000000000,14.0000000000,0.0020933151 +230.0000000000,15.0000000000,0.0021674652 +230.0000000000,16.0000000000,0.0017111567 +230.0000000000,17.0000000000,0.0012605521 +230.0000000000,18.0000000000,0.0011122519 +230.0000000000,19.0000000000,0.0008270591 +230.0000000000,20.0000000000,0.0007072781 +230.0000000000,21.0000000000,0.0005703856 +230.0000000000,22.0000000000,0.0003707506 +230.0000000000,23.0000000000,0.0001768195 +230.0000000000,24.0000000000,0.0001540041 +230.0000000000,25.0000000000,0.0000684463 +235.0000000000,0.0000000000,0.0000171116 +235.0000000000,1.0000000000,0.0002395619 +235.0000000000,2.0000000000,0.0004391969 +235.0000000000,3.0000000000,0.0006901666 +235.0000000000,4.0000000000,0.0007643167 +235.0000000000,5.0000000000,0.0011749943 +235.0000000000,6.0000000000,0.0014145562 +235.0000000000,7.0000000000,0.0016084873 +235.0000000000,8.0000000000,0.0019107917 +235.0000000000,9.0000000000,0.0018594570 +235.0000000000,10.0000000000,0.0020990189 +235.0000000000,11.0000000000,0.0022587269 +235.0000000000,12.0000000000,0.0020077572 +235.0000000000,13.0000000000,0.0018309377 +235.0000000000,14.0000000000,0.0019507187 +235.0000000000,15.0000000000,0.0013004791 +235.0000000000,16.0000000000,0.0013347023 +235.0000000000,17.0000000000,0.0011521789 +235.0000000000,18.0000000000,0.0009696555 +235.0000000000,19.0000000000,0.0008783938 +235.0000000000,20.0000000000,0.0005532740 +235.0000000000,21.0000000000,0.0003764545 +235.0000000000,22.0000000000,0.0002281542 +235.0000000000,23.0000000000,0.0001939311 +235.0000000000,24.0000000000,0.0000570386 +235.0000000000,25.0000000000,0.0000342231 +240.0000000000,0.0000000000,0.0000171116 +240.0000000000,1.0000000000,0.0001882272 +240.0000000000,2.0000000000,0.0003650468 +240.0000000000,3.0000000000,0.0006388319 +240.0000000000,4.0000000000,0.0009411362 +240.0000000000,5.0000000000,0.0009696555 +240.0000000000,6.0000000000,0.0012947753 +240.0000000000,7.0000000000,0.0014887064 +240.0000000000,8.0000000000,0.0016484143 +240.0000000000,9.0000000000,0.0020819074 +240.0000000000,10.0000000000,0.0021503536 +240.0000000000,11.0000000000,0.0018993840 +240.0000000000,12.0000000000,0.0020590919 +240.0000000000,13.0000000000,0.0018024184 +240.0000000000,14.0000000000,0.0017510837 +240.0000000000,15.0000000000,0.0012719598 +240.0000000000,16.0000000000,0.0011350673 +240.0000000000,17.0000000000,0.0009411362 +240.0000000000,18.0000000000,0.0009069131 +240.0000000000,19.0000000000,0.0006217203 +240.0000000000,20.0000000000,0.0004277892 +240.0000000000,21.0000000000,0.0003992699 +240.0000000000,22.0000000000,0.0001768195 +240.0000000000,23.0000000000,0.0001368925 +240.0000000000,24.0000000000,0.0000798540 +240.0000000000,25.0000000000,0.0000570386 +245.0000000000,0.0000000000,0.0000342231 +245.0000000000,1.0000000000,0.0002680812 +245.0000000000,2.0000000000,0.0004620123 +245.0000000000,3.0000000000,0.0007243897 +245.0000000000,4.0000000000,0.0008213552 +245.0000000000,5.0000000000,0.0009753593 +245.0000000000,6.0000000000,0.0011065480 +245.0000000000,7.0000000000,0.0013232945 +245.0000000000,8.0000000000,0.0015913758 +245.0000000000,9.0000000000,0.0019336071 +245.0000000000,10.0000000000,0.0018651608 +245.0000000000,11.0000000000,0.0018252339 +245.0000000000,12.0000000000,0.0016198950 +245.0000000000,13.0000000000,0.0016255989 +245.0000000000,14.0000000000,0.0014031485 +245.0000000000,15.0000000000,0.0012206251 +245.0000000000,16.0000000000,0.0010552133 +245.0000000000,17.0000000000,0.0009354324 +245.0000000000,18.0000000000,0.0009012092 +245.0000000000,19.0000000000,0.0006046087 +245.0000000000,20.0000000000,0.0004620123 +245.0000000000,21.0000000000,0.0003137121 +245.0000000000,22.0000000000,0.0001540041 +245.0000000000,23.0000000000,0.0000456308 +245.0000000000,24.0000000000,0.0000456308 +245.0000000000,25.0000000000,0.0000627424 +250.0000000000,0.0000000000,0.0000342231 +250.0000000000,1.0000000000,0.0002053388 +250.0000000000,2.0000000000,0.0004449008 +250.0000000000,3.0000000000,0.0007357974 +250.0000000000,4.0000000000,0.0008327629 +250.0000000000,5.0000000000,0.0011464750 +250.0000000000,6.0000000000,0.0012320329 +250.0000000000,7.0000000000,0.0015286334 +250.0000000000,8.0000000000,0.0016940452 +250.0000000000,9.0000000000,0.0017396760 +250.0000000000,10.0000000000,0.0018309377 +250.0000000000,11.0000000000,0.0017624914 +250.0000000000,12.0000000000,0.0014430755 +250.0000000000,13.0000000000,0.0014487794 +250.0000000000,14.0000000000,0.0012719598 +250.0000000000,15.0000000000,0.0010723249 +250.0000000000,16.0000000000,0.0009069131 +250.0000000000,17.0000000000,0.0007871321 +250.0000000000,18.0000000000,0.0006160164 +250.0000000000,19.0000000000,0.0004905316 +250.0000000000,20.0000000000,0.0004677162 +250.0000000000,21.0000000000,0.0002623774 +250.0000000000,22.0000000000,0.0002110427 +250.0000000000,23.0000000000,0.0000684463 +250.0000000000,24.0000000000,0.0000342231 +250.0000000000,25.0000000000,0.0000342231 +255.0000000000,0.0000000000,0.0000228154 +255.0000000000,1.0000000000,0.0002338581 +255.0000000000,2.0000000000,0.0004049738 +255.0000000000,3.0000000000,0.0006103126 +255.0000000000,4.0000000000,0.0008783938 +255.0000000000,5.0000000000,0.0009924709 +255.0000000000,6.0000000000,0.0012491444 +255.0000000000,7.0000000000,0.0013347023 +255.0000000000,8.0000000000,0.0015229295 +255.0000000000,9.0000000000,0.0015799681 +255.0000000000,10.0000000000,0.0016940452 +255.0000000000,11.0000000000,0.0013917408 +255.0000000000,12.0000000000,0.0015172256 +255.0000000000,13.0000000000,0.0013632215 +255.0000000000,14.0000000000,0.0009867671 +255.0000000000,15.0000000000,0.0010495095 +255.0000000000,16.0000000000,0.0009183208 +255.0000000000,17.0000000000,0.0008441707 +255.0000000000,18.0000000000,0.0005760894 +255.0000000000,19.0000000000,0.0004334930 +255.0000000000,20.0000000000,0.0003365275 +255.0000000000,21.0000000000,0.0001254848 +255.0000000000,22.0000000000,0.0001425964 +255.0000000000,23.0000000000,0.0000627424 +255.0000000000,24.0000000000,0.0000513347 +255.0000000000,25.0000000000,0.0000285193 +260.0000000000,0.0000000000,0.0000513347 +260.0000000000,1.0000000000,0.0002737851 +260.0000000000,2.0000000000,0.0004049738 +260.0000000000,3.0000000000,0.0007985398 +260.0000000000,4.0000000000,0.0008384668 +260.0000000000,5.0000000000,0.0011008442 +260.0000000000,6.0000000000,0.0012377367 +260.0000000000,7.0000000000,0.0013404061 +260.0000000000,8.0000000000,0.0012149213 +260.0000000000,9.0000000000,0.0015514488 +260.0000000000,10.0000000000,0.0013917408 +260.0000000000,11.0000000000,0.0012719598 +260.0000000000,12.0000000000,0.0013004791 +260.0000000000,13.0000000000,0.0011578827 +260.0000000000,14.0000000000,0.0012434406 +260.0000000000,15.0000000000,0.0010095825 +260.0000000000,16.0000000000,0.0006445357 +260.0000000000,17.0000000000,0.0005589779 +260.0000000000,18.0000000000,0.0004905316 +260.0000000000,19.0000000000,0.0003422313 +260.0000000000,20.0000000000,0.0002338581 +260.0000000000,21.0000000000,0.0001654118 +260.0000000000,22.0000000000,0.0001026694 +260.0000000000,23.0000000000,0.0000456308 +260.0000000000,24.0000000000,0.0000171116 +260.0000000000,25.0000000000,0.0000228154 +265.0000000000,0.0000000000,0.0000513347 +265.0000000000,1.0000000000,0.0001825234 +265.0000000000,2.0000000000,0.0003707506 +265.0000000000,3.0000000000,0.0007700205 +265.0000000000,4.0000000000,0.0008783938 +265.0000000000,5.0000000000,0.0009639516 +265.0000000000,6.0000000000,0.0012434406 +265.0000000000,7.0000000000,0.0013175907 +265.0000000000,8.0000000000,0.0013289984 +265.0000000000,9.0000000000,0.0015001141 +265.0000000000,10.0000000000,0.0012719598 +265.0000000000,11.0000000000,0.0013461100 +265.0000000000,12.0000000000,0.0012206251 +265.0000000000,13.0000000000,0.0010552133 +265.0000000000,14.0000000000,0.0009753593 +265.0000000000,15.0000000000,0.0007985398 +265.0000000000,16.0000000000,0.0006673511 +265.0000000000,17.0000000000,0.0005304586 +265.0000000000,18.0000000000,0.0004220853 +265.0000000000,19.0000000000,0.0004506046 +265.0000000000,20.0000000000,0.0002452658 +265.0000000000,21.0000000000,0.0001711157 +265.0000000000,22.0000000000,0.0001254848 +265.0000000000,23.0000000000,0.0000456308 +265.0000000000,24.0000000000,0.0000114077 +265.0000000000,25.0000000000,0.0000057039 +270.0000000000,0.0000000000,0.0000171116 +270.0000000000,1.0000000000,0.0001825234 +270.0000000000,2.0000000000,0.0004220853 +270.0000000000,3.0000000000,0.0006730550 +270.0000000000,4.0000000000,0.0007928360 +270.0000000000,5.0000000000,0.0011350673 +270.0000000000,6.0000000000,0.0011464750 +270.0000000000,7.0000000000,0.0013175907 +270.0000000000,8.0000000000,0.0012662560 +270.0000000000,9.0000000000,0.0014145562 +270.0000000000,10.0000000000,0.0014088524 +270.0000000000,11.0000000000,0.0012263290 +270.0000000000,12.0000000000,0.0010552133 +270.0000000000,13.0000000000,0.0011008442 +270.0000000000,14.0000000000,0.0009468401 +270.0000000000,15.0000000000,0.0008555784 +270.0000000000,16.0000000000,0.0006958704 +270.0000000000,17.0000000000,0.0006046087 +270.0000000000,18.0000000000,0.0004391969 +270.0000000000,19.0000000000,0.0003194159 +270.0000000000,20.0000000000,0.0002851928 +270.0000000000,21.0000000000,0.0002110427 +270.0000000000,22.0000000000,0.0001197810 +270.0000000000,23.0000000000,0.0000285193 +270.0000000000,24.0000000000,0.0000171116 +270.0000000000,25.0000000000,0.0000171116 +275.0000000000,0.0000000000,0.0000513347 +275.0000000000,1.0000000000,0.0002509697 +275.0000000000,2.0000000000,0.0005133470 +275.0000000000,3.0000000000,0.0005646817 +275.0000000000,4.0000000000,0.0008898015 +275.0000000000,5.0000000000,0.0010495095 +275.0000000000,6.0000000000,0.0012035136 +275.0000000000,7.0000000000,0.0013061830 +275.0000000000,8.0000000000,0.0014772987 +275.0000000000,9.0000000000,0.0015058179 +275.0000000000,10.0000000000,0.0015514488 +275.0000000000,11.0000000000,0.0013118868 +275.0000000000,12.0000000000,0.0010951403 +275.0000000000,13.0000000000,0.0010438056 +275.0000000000,14.0000000000,0.0007529090 +275.0000000000,15.0000000000,0.0007243897 +275.0000000000,16.0000000000,0.0006445357 +275.0000000000,17.0000000000,0.0004734200 +275.0000000000,18.0000000000,0.0004620123 +275.0000000000,19.0000000000,0.0003707506 +275.0000000000,20.0000000000,0.0001996350 +275.0000000000,21.0000000000,0.0001425964 +275.0000000000,22.0000000000,0.0000912617 +275.0000000000,23.0000000000,0.0000285193 +275.0000000000,24.0000000000,0.0000399270 +275.0000000000,25.0000000000,0.0000171116 +280.0000000000,0.0000000000,0.0000513347 +280.0000000000,1.0000000000,0.0002167465 +280.0000000000,2.0000000000,0.0003194159 +280.0000000000,3.0000000000,0.0006616473 +280.0000000000,4.0000000000,0.0007472051 +280.0000000000,5.0000000000,0.0009069131 +280.0000000000,6.0000000000,0.0010837326 +280.0000000000,7.0000000000,0.0015343372 +280.0000000000,8.0000000000,0.0014430755 +280.0000000000,9.0000000000,0.0015571526 +280.0000000000,10.0000000000,0.0013860370 +280.0000000000,11.0000000000,0.0013860370 +280.0000000000,12.0000000000,0.0011464750 +280.0000000000,13.0000000000,0.0010095825 +280.0000000000,14.0000000000,0.0007814282 +280.0000000000,15.0000000000,0.0007700205 +280.0000000000,16.0000000000,0.0005475702 +280.0000000000,17.0000000000,0.0005361624 +280.0000000000,18.0000000000,0.0004049738 +280.0000000000,19.0000000000,0.0003422313 +280.0000000000,20.0000000000,0.0002395619 +280.0000000000,21.0000000000,0.0001825234 +280.0000000000,22.0000000000,0.0000798540 +280.0000000000,23.0000000000,0.0000399270 +280.0000000000,24.0000000000,0.0000399270 +280.0000000000,25.0000000000,0.0000171116 +285.0000000000,0.0000000000,0.0000228154 +285.0000000000,1.0000000000,0.0002224504 +285.0000000000,2.0000000000,0.0003764545 +285.0000000000,3.0000000000,0.0005304586 +285.0000000000,4.0000000000,0.0007300935 +285.0000000000,5.0000000000,0.0010266940 +285.0000000000,6.0000000000,0.0012092174 +285.0000000000,7.0000000000,0.0013175907 +285.0000000000,8.0000000000,0.0013461100 +285.0000000000,9.0000000000,0.0015970796 +285.0000000000,10.0000000000,0.0012947753 +285.0000000000,11.0000000000,0.0012035136 +285.0000000000,12.0000000000,0.0011464750 +285.0000000000,13.0000000000,0.0007928360 +285.0000000000,14.0000000000,0.0008270591 +285.0000000000,15.0000000000,0.0006901666 +285.0000000000,16.0000000000,0.0005304586 +285.0000000000,17.0000000000,0.0003878622 +285.0000000000,18.0000000000,0.0003479352 +285.0000000000,19.0000000000,0.0002966005 +285.0000000000,20.0000000000,0.0002110427 +285.0000000000,21.0000000000,0.0001540041 +285.0000000000,22.0000000000,0.0000741501 +285.0000000000,23.0000000000,0.0000399270 +285.0000000000,24.0000000000,0.0000285193 +285.0000000000,25.0000000000,0.0000000000 +290.0000000000,0.0000000000,0.0000342231 +290.0000000000,1.0000000000,0.0002395619 +290.0000000000,2.0000000000,0.0004563085 +290.0000000000,3.0000000000,0.0005190509 +290.0000000000,4.0000000000,0.0007928360 +290.0000000000,5.0000000000,0.0009525439 +290.0000000000,6.0000000000,0.0010837326 +290.0000000000,7.0000000000,0.0013347023 +290.0000000000,8.0000000000,0.0014430755 +290.0000000000,9.0000000000,0.0012320329 +290.0000000000,10.0000000000,0.0013746292 +290.0000000000,11.0000000000,0.0010780287 +290.0000000000,12.0000000000,0.0009354324 +290.0000000000,13.0000000000,0.0009582478 +290.0000000000,14.0000000000,0.0006958704 +290.0000000000,15.0000000000,0.0005703856 +290.0000000000,16.0000000000,0.0003764545 +290.0000000000,17.0000000000,0.0003878622 +290.0000000000,18.0000000000,0.0002167465 +290.0000000000,19.0000000000,0.0001882272 +290.0000000000,20.0000000000,0.0001711157 +290.0000000000,21.0000000000,0.0000969655 +290.0000000000,22.0000000000,0.0000684463 +290.0000000000,23.0000000000,0.0000228154 +290.0000000000,24.0000000000,0.0000114077 +290.0000000000,25.0000000000,0.0000171116 +295.0000000000,0.0000000000,0.0000399270 +295.0000000000,1.0000000000,0.0002281542 +295.0000000000,2.0000000000,0.0004220853 +295.0000000000,3.0000000000,0.0006388319 +295.0000000000,4.0000000000,0.0007757244 +295.0000000000,5.0000000000,0.0008955054 +295.0000000000,6.0000000000,0.0012092174 +295.0000000000,7.0000000000,0.0013232945 +295.0000000000,8.0000000000,0.0011978097 +295.0000000000,9.0000000000,0.0013860370 +295.0000000000,10.0000000000,0.0013461100 +295.0000000000,11.0000000000,0.0013404061 +295.0000000000,12.0000000000,0.0010552133 +295.0000000000,13.0000000000,0.0009525439 +295.0000000000,14.0000000000,0.0006445357 +295.0000000000,15.0000000000,0.0005361624 +295.0000000000,16.0000000000,0.0002966005 +295.0000000000,17.0000000000,0.0002851928 +295.0000000000,18.0000000000,0.0001825234 +295.0000000000,19.0000000000,0.0001483003 +295.0000000000,20.0000000000,0.0000969655 +295.0000000000,21.0000000000,0.0000570386 +295.0000000000,22.0000000000,0.0000342231 +295.0000000000,23.0000000000,0.0000171116 +295.0000000000,24.0000000000,0.0000114077 +295.0000000000,25.0000000000,0.0000114077 +300.0000000000,0.0000000000,0.0000684463 +300.0000000000,1.0000000000,0.0002281542 +300.0000000000,2.0000000000,0.0004220853 +300.0000000000,3.0000000000,0.0006046087 +300.0000000000,4.0000000000,0.0008327629 +300.0000000000,5.0000000000,0.0009126169 +300.0000000000,6.0000000000,0.0011978097 +300.0000000000,7.0000000000,0.0011978097 +300.0000000000,8.0000000000,0.0013746292 +300.0000000000,9.0000000000,0.0013974447 +300.0000000000,10.0000000000,0.0012206251 +300.0000000000,11.0000000000,0.0011635866 +300.0000000000,12.0000000000,0.0009753593 +300.0000000000,13.0000000000,0.0007757244 +300.0000000000,14.0000000000,0.0005475702 +300.0000000000,15.0000000000,0.0004563085 +300.0000000000,16.0000000000,0.0003935661 +300.0000000000,17.0000000000,0.0002794889 +300.0000000000,18.0000000000,0.0001996350 +300.0000000000,19.0000000000,0.0001368925 +300.0000000000,20.0000000000,0.0000513347 +300.0000000000,21.0000000000,0.0000627424 +300.0000000000,22.0000000000,0.0000285193 +300.0000000000,23.0000000000,0.0000228154 +300.0000000000,24.0000000000,0.0000171116 +300.0000000000,25.0000000000,0.0000057039 +305.0000000000,0.0000000000,0.0000285193 +305.0000000000,1.0000000000,0.0001996350 +305.0000000000,2.0000000000,0.0004163815 +305.0000000000,3.0000000000,0.0006217203 +305.0000000000,4.0000000000,0.0007757244 +305.0000000000,5.0000000000,0.0008783938 +305.0000000000,6.0000000000,0.0011350673 +305.0000000000,7.0000000000,0.0012890714 +305.0000000000,8.0000000000,0.0012491444 +305.0000000000,9.0000000000,0.0013917408 +305.0000000000,10.0000000000,0.0012662560 +305.0000000000,11.0000000000,0.0011293634 +305.0000000000,12.0000000000,0.0008898015 +305.0000000000,13.0000000000,0.0008612822 +305.0000000000,14.0000000000,0.0006274241 +305.0000000000,15.0000000000,0.0005190509 +305.0000000000,16.0000000000,0.0002737851 +305.0000000000,17.0000000000,0.0002509697 +305.0000000000,18.0000000000,0.0001597080 +305.0000000000,19.0000000000,0.0000627424 +305.0000000000,20.0000000000,0.0001026694 +305.0000000000,21.0000000000,0.0000513347 +305.0000000000,22.0000000000,0.0000570386 +305.0000000000,23.0000000000,0.0000057039 +305.0000000000,24.0000000000,0.0000000000 +305.0000000000,25.0000000000,0.0000000000 +310.0000000000,0.0000000000,0.0000342231 +310.0000000000,1.0000000000,0.0002167465 +310.0000000000,2.0000000000,0.0004563085 +310.0000000000,3.0000000000,0.0006103126 +310.0000000000,4.0000000000,0.0008327629 +310.0000000000,5.0000000000,0.0009753593 +310.0000000000,6.0000000000,0.0010209902 +310.0000000000,7.0000000000,0.0013803331 +310.0000000000,8.0000000000,0.0012662560 +310.0000000000,9.0000000000,0.0012491444 +310.0000000000,10.0000000000,0.0012605521 +310.0000000000,11.0000000000,0.0009924709 +310.0000000000,12.0000000000,0.0009981748 +310.0000000000,13.0000000000,0.0008156514 +310.0000000000,14.0000000000,0.0006901666 +310.0000000000,15.0000000000,0.0004848277 +310.0000000000,16.0000000000,0.0004677162 +310.0000000000,17.0000000000,0.0002395619 +310.0000000000,18.0000000000,0.0002224504 +310.0000000000,19.0000000000,0.0001311887 +310.0000000000,20.0000000000,0.0000855578 +310.0000000000,21.0000000000,0.0000798540 +310.0000000000,22.0000000000,0.0000741501 +310.0000000000,23.0000000000,0.0000285193 +310.0000000000,24.0000000000,0.0000171116 +310.0000000000,25.0000000000,0.0000000000 +315.0000000000,0.0000000000,0.0000285193 +315.0000000000,1.0000000000,0.0001939311 +315.0000000000,2.0000000000,0.0003992699 +315.0000000000,3.0000000000,0.0005532740 +315.0000000000,4.0000000000,0.0008441707 +315.0000000000,5.0000000000,0.0011008442 +315.0000000000,6.0000000000,0.0012263290 +315.0000000000,7.0000000000,0.0012491444 +315.0000000000,8.0000000000,0.0012263290 +315.0000000000,9.0000000000,0.0012320329 +315.0000000000,10.0000000000,0.0011978097 +315.0000000000,11.0000000000,0.0008213552 +315.0000000000,12.0000000000,0.0008955054 +315.0000000000,13.0000000000,0.0008213552 +315.0000000000,14.0000000000,0.0007186858 +315.0000000000,15.0000000000,0.0004848277 +315.0000000000,16.0000000000,0.0003308236 +315.0000000000,17.0000000000,0.0002966005 +315.0000000000,18.0000000000,0.0002224504 +315.0000000000,19.0000000000,0.0001425964 +315.0000000000,20.0000000000,0.0001254848 +315.0000000000,21.0000000000,0.0000912617 +315.0000000000,22.0000000000,0.0000456308 +315.0000000000,23.0000000000,0.0000399270 +315.0000000000,24.0000000000,0.0000114077 +315.0000000000,25.0000000000,0.0000000000 +320.0000000000,0.0000000000,0.0000741501 +320.0000000000,1.0000000000,0.0002509697 +320.0000000000,2.0000000000,0.0003878622 +320.0000000000,3.0000000000,0.0006559434 +320.0000000000,4.0000000000,0.0007015743 +320.0000000000,5.0000000000,0.0009582478 +320.0000000000,6.0000000000,0.0009810632 +320.0000000000,7.0000000000,0.0011806982 +320.0000000000,8.0000000000,0.0010666210 +320.0000000000,9.0000000000,0.0010152863 +320.0000000000,10.0000000000,0.0012206251 +320.0000000000,11.0000000000,0.0010152863 +320.0000000000,12.0000000000,0.0008327629 +320.0000000000,13.0000000000,0.0007129820 +320.0000000000,14.0000000000,0.0005760894 +320.0000000000,15.0000000000,0.0005532740 +320.0000000000,16.0000000000,0.0003593429 +320.0000000000,17.0000000000,0.0003137121 +320.0000000000,18.0000000000,0.0001483003 +320.0000000000,19.0000000000,0.0001425964 +320.0000000000,20.0000000000,0.0000912617 +320.0000000000,21.0000000000,0.0000456308 +320.0000000000,22.0000000000,0.0000342231 +320.0000000000,23.0000000000,0.0000228154 +320.0000000000,24.0000000000,0.0000171116 +320.0000000000,25.0000000000,0.0000000000 +325.0000000000,0.0000000000,0.0000570386 +325.0000000000,1.0000000000,0.0002110427 +325.0000000000,2.0000000000,0.0004391969 +325.0000000000,3.0000000000,0.0006160164 +325.0000000000,4.0000000000,0.0008270591 +325.0000000000,5.0000000000,0.0008955054 +325.0000000000,6.0000000000,0.0011293634 +325.0000000000,7.0000000000,0.0011236596 +325.0000000000,8.0000000000,0.0011464750 +325.0000000000,9.0000000000,0.0011179557 +325.0000000000,10.0000000000,0.0010837326 +325.0000000000,11.0000000000,0.0009639516 +325.0000000000,12.0000000000,0.0009753593 +325.0000000000,13.0000000000,0.0006958704 +325.0000000000,14.0000000000,0.0006331280 +325.0000000000,15.0000000000,0.0004734200 +325.0000000000,16.0000000000,0.0003479352 +325.0000000000,17.0000000000,0.0002281542 +325.0000000000,18.0000000000,0.0002053388 +325.0000000000,19.0000000000,0.0002110427 +325.0000000000,20.0000000000,0.0001140771 +325.0000000000,21.0000000000,0.0000627424 +325.0000000000,22.0000000000,0.0000456308 +325.0000000000,23.0000000000,0.0000114077 +325.0000000000,24.0000000000,0.0000057039 +325.0000000000,25.0000000000,0.0000000000 +330.0000000000,0.0000000000,0.0000456308 +330.0000000000,1.0000000000,0.0002167465 +330.0000000000,2.0000000000,0.0004848277 +330.0000000000,3.0000000000,0.0006673511 +330.0000000000,4.0000000000,0.0007757244 +330.0000000000,5.0000000000,0.0009525439 +330.0000000000,6.0000000000,0.0010837326 +330.0000000000,7.0000000000,0.0011978097 +330.0000000000,8.0000000000,0.0009525439 +330.0000000000,9.0000000000,0.0009240246 +330.0000000000,10.0000000000,0.0008726899 +330.0000000000,11.0000000000,0.0008498745 +330.0000000000,12.0000000000,0.0008669861 +330.0000000000,13.0000000000,0.0006901666 +330.0000000000,14.0000000000,0.0005760894 +330.0000000000,15.0000000000,0.0004620123 +330.0000000000,16.0000000000,0.0003422313 +330.0000000000,17.0000000000,0.0002566735 +330.0000000000,18.0000000000,0.0002623774 +330.0000000000,19.0000000000,0.0001254848 +330.0000000000,20.0000000000,0.0001254848 +330.0000000000,21.0000000000,0.0000399270 +330.0000000000,22.0000000000,0.0000342231 +330.0000000000,23.0000000000,0.0000000000 +330.0000000000,24.0000000000,0.0000000000 +330.0000000000,25.0000000000,0.0000000000 +335.0000000000,0.0000000000,0.0000342231 +335.0000000000,1.0000000000,0.0001654118 +335.0000000000,2.0000000000,0.0005076432 +335.0000000000,3.0000000000,0.0006673511 +335.0000000000,4.0000000000,0.0007814282 +335.0000000000,5.0000000000,0.0010666210 +335.0000000000,6.0000000000,0.0009696555 +335.0000000000,7.0000000000,0.0012320329 +335.0000000000,8.0000000000,0.0011749943 +335.0000000000,9.0000000000,0.0009924709 +335.0000000000,10.0000000000,0.0010323979 +335.0000000000,11.0000000000,0.0006901666 +335.0000000000,12.0000000000,0.0007529090 +335.0000000000,13.0000000000,0.0005589779 +335.0000000000,14.0000000000,0.0004220853 +335.0000000000,15.0000000000,0.0003365275 +335.0000000000,16.0000000000,0.0002452658 +335.0000000000,17.0000000000,0.0002053388 +335.0000000000,18.0000000000,0.0001140771 +335.0000000000,19.0000000000,0.0000912617 +335.0000000000,20.0000000000,0.0000684463 +335.0000000000,21.0000000000,0.0000171116 +335.0000000000,22.0000000000,0.0000114077 +335.0000000000,23.0000000000,0.0000000000 +335.0000000000,24.0000000000,0.0000000000 +335.0000000000,25.0000000000,0.0000000000 +340.0000000000,0.0000000000,0.0000285193 +340.0000000000,1.0000000000,0.0002395619 +340.0000000000,2.0000000000,0.0003536391 +340.0000000000,3.0000000000,0.0005019393 +340.0000000000,4.0000000000,0.0007700205 +340.0000000000,5.0000000000,0.0008669861 +340.0000000000,6.0000000000,0.0008555784 +340.0000000000,7.0000000000,0.0010837326 +340.0000000000,8.0000000000,0.0010095825 +340.0000000000,9.0000000000,0.0010266940 +340.0000000000,10.0000000000,0.0010266940 +340.0000000000,11.0000000000,0.0009012092 +340.0000000000,12.0000000000,0.0008270591 +340.0000000000,13.0000000000,0.0006559434 +340.0000000000,14.0000000000,0.0004449008 +340.0000000000,15.0000000000,0.0003365275 +340.0000000000,16.0000000000,0.0002053388 +340.0000000000,17.0000000000,0.0001425964 +340.0000000000,18.0000000000,0.0000855578 +340.0000000000,19.0000000000,0.0000855578 +340.0000000000,20.0000000000,0.0000342231 +340.0000000000,21.0000000000,0.0000057039 +340.0000000000,22.0000000000,0.0000057039 +340.0000000000,23.0000000000,0.0000000000 +340.0000000000,24.0000000000,0.0000000000 +340.0000000000,25.0000000000,0.0000000000 +345.0000000000,0.0000000000,0.0000342231 +345.0000000000,1.0000000000,0.0001825234 +345.0000000000,2.0000000000,0.0003764545 +345.0000000000,3.0000000000,0.0006103126 +345.0000000000,4.0000000000,0.0007472051 +345.0000000000,5.0000000000,0.0008726899 +345.0000000000,6.0000000000,0.0008726899 +345.0000000000,7.0000000000,0.0009297285 +345.0000000000,8.0000000000,0.0010381018 +345.0000000000,9.0000000000,0.0008099475 +345.0000000000,10.0000000000,0.0009753593 +345.0000000000,11.0000000000,0.0008612822 +345.0000000000,12.0000000000,0.0005989049 +345.0000000000,13.0000000000,0.0005703856 +345.0000000000,14.0000000000,0.0003308236 +345.0000000000,15.0000000000,0.0003194159 +345.0000000000,16.0000000000,0.0002281542 +345.0000000000,17.0000000000,0.0000798540 +345.0000000000,18.0000000000,0.0000855578 +345.0000000000,19.0000000000,0.0000399270 +345.0000000000,20.0000000000,0.0000228154 +345.0000000000,21.0000000000,0.0000000000 +345.0000000000,22.0000000000,0.0000000000 +345.0000000000,23.0000000000,0.0000000000 +345.0000000000,24.0000000000,0.0000000000 +345.0000000000,25.0000000000,0.0000000000 +350.0000000000,0.0000000000,0.0000285193 +350.0000000000,1.0000000000,0.0001711157 +350.0000000000,2.0000000000,0.0004506046 +350.0000000000,3.0000000000,0.0005532740 +350.0000000000,4.0000000000,0.0006844627 +350.0000000000,5.0000000000,0.0007015743 +350.0000000000,6.0000000000,0.0008270591 +350.0000000000,7.0000000000,0.0008840977 +350.0000000000,8.0000000000,0.0009981748 +350.0000000000,9.0000000000,0.0009297285 +350.0000000000,10.0000000000,0.0009582478 +350.0000000000,11.0000000000,0.0007015743 +350.0000000000,12.0000000000,0.0005703856 +350.0000000000,13.0000000000,0.0005817933 +350.0000000000,14.0000000000,0.0002110427 +350.0000000000,15.0000000000,0.0001654118 +350.0000000000,16.0000000000,0.0001711157 +350.0000000000,17.0000000000,0.0001026694 +350.0000000000,18.0000000000,0.0000627424 +350.0000000000,19.0000000000,0.0000228154 +350.0000000000,20.0000000000,0.0000114077 +350.0000000000,21.0000000000,0.0000114077 +350.0000000000,22.0000000000,0.0000000000 +350.0000000000,23.0000000000,0.0000000000 +350.0000000000,24.0000000000,0.0000000000 +350.0000000000,25.0000000000,0.0000000000 +355.0000000000,0.0000000000,0.0000342231 +355.0000000000,1.0000000000,0.0002167465 +355.0000000000,2.0000000000,0.0003650468 +355.0000000000,3.0000000000,0.0005418663 +355.0000000000,4.0000000000,0.0006160164 +355.0000000000,5.0000000000,0.0008955054 +355.0000000000,6.0000000000,0.0008498745 +355.0000000000,7.0000000000,0.0009126169 +355.0000000000,8.0000000000,0.0010438056 +355.0000000000,9.0000000000,0.0008898015 +355.0000000000,10.0000000000,0.0009240246 +355.0000000000,11.0000000000,0.0007415013 +355.0000000000,12.0000000000,0.0005874971 +355.0000000000,13.0000000000,0.0003365275 +355.0000000000,14.0000000000,0.0003137121 +355.0000000000,15.0000000000,0.0002338581 +355.0000000000,16.0000000000,0.0001425964 +355.0000000000,17.0000000000,0.0001140771 +355.0000000000,18.0000000000,0.0000627424 +355.0000000000,19.0000000000,0.0000114077 +355.0000000000,20.0000000000,0.0000171116 +355.0000000000,21.0000000000,0.0000000000 +355.0000000000,22.0000000000,0.0000057039 +355.0000000000,23.0000000000,0.0000000000 +355.0000000000,24.0000000000,0.0000000000 +355.0000000000,25.0000000000,0.0000000000 diff --git a/flowers/__init__.py b/flowers/__init__.py new file mode 100644 index 00000000..03bb1b92 --- /dev/null +++ b/flowers/__init__.py @@ -0,0 +1 @@ +from .flowers_model import FlowersModel diff --git a/flowers/flowers_model.py b/flowers/flowers_model.py new file mode 100644 index 00000000..349c6880 --- /dev/null +++ b/flowers/flowers_model.py @@ -0,0 +1,369 @@ +import copy + +import numpy as np +import pandas as pd + +import flowers.tools as tl + + +class FlowersModel: + """ + FlowersModel is a high-level user interface to the FLOWERS AEP model. + + Args: + wind_rose (pandas.DataFrame): A dataframe for the wind rose in the FLORIS + format containing the following information: + - 'ws' (float): wind speeds [m/s] + - 'wd' (float): wind directions [deg] + - 'freq_val' (float): frequency for each wind speed and direction + layout_x (numpy.array(float)): x-positions of each turbine [m] + layout_y (numpy.array(float)): y-positions of each turbine [m] + num_terms (int, optional): number of Fourier modes + k (float, optional): wake expansion rate + turbine (str, optional): turbine type: + - 'nrel_5MW' (default) + + """ + + ########################################################################### + # Initialization tools + ########################################################################### + + def __init__( + self, wind_rose, layout_x, layout_y, num_terms=0, k=0.05, turbine=None + ): + + self.wind_rose = wind_rose + self.layout_x = layout_x + self.layout_y = layout_y + self.k = k + + if turbine is None: + turbine = "nrel_5MW" + if type(turbine) == str: + if turbine == "nrel_5MW": + self.turbine = "nrel_5MW" + self.D = 126.0 + self.U = 25.0 + elif turbine == "iea_22MW": + self.turbine = "iea_22MW" + self.D = 283.2 + self.U = 25.0 + else: + raise NotImplementedError( + f"turbine type ({turbine}) has not been implemented" + ) + + self._fourier_coefficients(num_terms=num_terms) + + elif type(turbine) == dict: + assert "ct" in turbine + assert "cp" in turbine + assert "D" in turbine + assert "U" in turbine + self.D = turbine["D"] + self.U = turbine["U"] + ct = turbine["ct"] + cp = turbine["cp"] + u_ct = turbine["u_ct"] + u_cp = turbine["u_cp"] + + self._fourier_coefficients( + num_terms=num_terms, + u_ct_table=u_ct, + u_cp_table=u_cp, + ct_table=ct, + cp_table=cp, + ) + + else: + raise IOError("supplied turbine type was not recognized!") + + def reinitialize( + self, wind_rose=None, layout_x=None, layout_y=None, num_terms=None, k=None + ): + + if wind_rose is not None: + self.wind_rose = wind_rose + self._fourier_coefficients(num_terms=num_terms) + + if num_terms is not None: + self._fourier_coefficients(num_terms=num_terms) + + if layout_x is not None: + self.layout_x = layout_x + + if layout_y is not None: + self.layout_y = layout_y + + if k is not None: + self.k = k + + ########################################################################### + # User functions + ########################################################################### + + def get_layout(self): + return self.layout_x, self.layout_y + + def get_wind_rose(self): + return self.wind_rose + + def get_num_modes(self): + return len(self.fs) + + def calculate_aep(self, rho_density=1.225, gradient=False): + """ + Compute farm AEP (and Cartesian gradients) for the given layout and wind rose. + + Returns: + aep (float): farm AEP [Wh] + gradient (numpy.array(float)): (dAEP/dx, dAEP/dy) for each turbine [Wh/m] + """ + + # Power component from freestream + u0 = self.fs.c[0] + + # Normalize and reshape relative positions into symmetric 2D array + xx = (self.layout_x - np.reshape(self.layout_x, (-1, 1))) / self.D + yy = (self.layout_y - np.reshape(self.layout_y, (-1, 1))) / self.D + + # Convert to normalized polar coordinates + R = np.sqrt(xx**2 + yy**2) + THETA = np.arctan2(yy, xx) / (2 * np.pi) + + # Set up mask for rotor swept area + mask_area = np.array(R <= 0.5, dtype=int) + mask_val = self.fs.c[0] + + # Critical polar angle of wake edge (as a function of distance from turbine) + with np.errstate(divide="ignore", invalid="ignore"): + theta_c = np.arctan( + (1 / (2 * R) + self.k * np.sqrt(1 + self.k**2 - (2 * R) ** (-2))) + / (-self.k / (2 * R) + np.sqrt(1 + self.k**2 - (2 * R) ** (-2))) + ) / (2 * np.pi) + theta_c = np.nan_to_num(theta_c) + + # Contribution from zero-frequency Fourier mode + du = ( + self.fs.a[0] + * theta_c + / (2 * self.k * R + 1) ** 2 + * ( + 1 + + (8 * np.pi**2 * theta_c**2 * self.k * R) / (3 * (2 * self.k * R + 1)) + ) + ) + + # Initialize gradient and calculate zero-frequency modes + if gradient == True: + grad = np.zeros((len(self.layout_x), 2)) + + # Change in theta_c wrt radius + dtdr = -1 / (4 * np.pi * R**2 * np.sqrt(self.k**2 - (2 * R) ** (-2) + 1)) + dtdr = np.nan_to_num(dtdr) + + # Zero-frequency mode of change in power deficit wrt radius + dpdr = ( + -4 + * self.fs.a[0] + * self.k + * theta_c + * ( + 3 + + 6 * self.k * R + + 2 * np.pi**2 * (4 * self.k * R - 1) * theta_c**2 + ) + + 3 + * self.fs.a[0] + * (1 + 2 * self.k * R) + * (1 + 2 * self.k * R + 8 * np.pi**2 * self.k * R * theta_c**2) + * dtdr + ) / (3 * (1 + 2 * self.k * R) ** 4) + + for m in np.arange(1, len(self.fs.b)): + du += ( + 1 + / (np.pi * m * (2 * self.k * R + 1) ** 2) + * ( + self.fs.a[m] * np.cos(2 * np.pi * m * THETA) + + self.fs.b[m] * np.sin(2 * np.pi * m * THETA) + ) + * ( + np.sin(2 * np.pi * m * theta_c) + + 2 + * self.k + * R + / (m**2 * (2 * self.k * R + 1)) + * ( + ((2 * np.pi * theta_c * m) ** 2 - 2) + * np.sin(2 * np.pi * m * theta_c) + + 4 * np.pi * m * theta_c * np.cos(2 * np.pi * m * theta_c) + ) + ) + ) + + if gradient == True: + dpdt = 0 + for m in np.arange(1, len(self.fs.b)): + # Higher Fourier modes of change in power deficit wrt angle + dpdt += ( + 2 + / (2 * self.k * R + 1) ** 2 + * ( + self.fs.b[m] * np.cos(2 * np.pi * m * THETA) + - self.fs.a[m] * np.sin(2 * np.pi * m * THETA) + ) + * ( + np.sin(2 * np.pi * m * theta_c) + + 2 + * self.k + * R + / (m**2 * (2 * self.k * R + 1)) + * ( + ((2 * np.pi * theta_c * m) ** 2 - 2) + * np.sin(2 * np.pi * m * theta_c) + + 4 * np.pi * m * theta_c * np.cos(2 * np.pi * m * theta_c) + ) + ) + ) + + # Higher Fourier modes of change in power deficit wrt radius + dpdr += ( + ( + self.fs.a[m] * np.cos(2 * np.pi * m * THETA) + + self.fs.b[m] * np.sin(2 * np.pi * m * THETA) + ) + / (np.pi * m**3 * (2 * self.k * R + 1) ** 4) + * ( + -4 + * self.k + * np.sin(2 * np.pi * m * theta_c) + * ( + 1 + + m**2 + + 2 * self.k * R * (m**2 - 2) + + 2 * np.pi**2 * m**2 * (4 * self.k * R - 1) * theta_c**2 + ) + + 2 + * np.pi + * m + * np.cos(2 * np.pi * m * theta_c) + * ( + 4 * self.k * (1 - 4 * self.k * R) * theta_c + + m**2 + * (2 * self.k * R + 1) + * ( + 1 + + 2 * self.k * R + + 8 * np.pi**2 * self.k * R * theta_c**2 + ) + * dtdr + ) + ) + ) + + # Apply mask for points within rotor radius + du = du * (1 - mask_area) + mask_val * mask_area + np.fill_diagonal(du, 0.0) + + # Sum power for each turbine + du = np.sum(du, axis=1) + aep = np.sum((u0 - du) ** 3) + aep *= np.pi / 8 * rho_density * self.D**2 * self.U**3 * 8760 + + # Complete gradient calculation + if gradient == True: + dx = xx / R * dpdr + -yy / (2 * np.pi * R**2) * dpdt + dy = yy / R * dpdr + xx / (2 * np.pi * R**2) * dpdt + + dx = np.nan_to_num(dx) + dy = np.nan_to_num(dy) + + coeff = (u0 - du) ** 2 + for i in range(len(grad)): + # Isolate gradient to turbine 'i' + grad_mask = np.zeros_like(xx) + grad_mask[i, :] = -1.0 + grad_mask[:, i] = 1.0 + + grad[i, 0] = np.sum(coeff * np.sum(dx * grad_mask, axis=1)) + grad[i, 1] = np.sum(coeff * np.sum(dy * grad_mask, axis=1)) + + grad *= -3 * np.pi / 8 * rho_density * self.D * self.U**3 * 8760 + + return aep, grad + + else: + return aep + + ########################################################################### + # Private functions + ########################################################################### + + def _fourier_coefficients( + self, + num_terms=0, + ct_table=None, + cp_table=None, + u_ct_table=None, + u_cp_table=None, + ): + """ + Compute the Fourier series expansion coefficients from the wind rose. + Modifies the FlowersModel in place to add a Fourier coefficients + dataframe: + fs (pandas:dataframe): Fourier coefficients used to expand the wind rose: + - 'a_free': real coefficients of freestream component + - 'a_wake': real coefficients of wake component + - 'b_wake': imaginary coefficients of wake component + + Args: + num_terms (int, optional): the number of Fourier modes to save in the range + [1, floor(num_wind_directions/2)] + + """ + + # Resample wind rose for average wind speed per wind direction + wr = copy.deepcopy(self.wind_rose) + wr = tl.resample_average_ws_by_wd(wr) + + # Transform wind direction to polar angle + wr["wd"] = np.remainder(450 - wr.wd, 360) + wr.sort_values("wd", inplace=True) + # wr.loc[len(wr)] = wr.iloc[0] + # wr.freq_val /= np.sum(wr.freq_val) + + # Normalize wind speed by cut-out speed + wr["ws"] /= self.U + u_ct_table = None if u_ct_table is None else np.array(u_ct_table) / self.U + u_cp_table = None if u_cp_table is None else np.array(u_cp_table) / self.U + + # Look up thrust and power coefficients for each wind direction bin + if (u_ct_table is None) or (ct_table is None): + ct = tl.ct_lookup(wr.ws, self.turbine) + else: + # pull out the thrust lookup interpolation + ct = np.interp(wr.ws, u_ct_table, ct_table) + if (u_cp_table is None) or (cp_table is None): + cp = tl.cp_lookup(wr.ws, self.turbine) + else: + # pull out the power lookup interpolation + cp = np.interp(wr.ws, u_cp_table, cp_table) + + # Average freestream term + c = np.sum(cp ** (1 / 3) * wr.ws * wr.freq_val) + + # Fourier expansion of wake deficit term + c1 = cp ** (1 / 3) * (1 - np.sqrt(1 - ct)) * wr.ws * wr.freq_val + c1ft = 2 * np.fft.rfft(c1) + a = c1ft.real + b = -c1ft.imag + + # Truncate Fourier series to specified number of modes + if num_terms > 0 and num_terms <= len(a): + a = a[0:num_terms] + b = b[0:num_terms] + + # Compile Fourier coefficients + self.fs = pd.DataFrame({"a": a, "b": b, "c": c}) diff --git a/flowers/optimization/__init__.py b/flowers/optimization/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/flowers/optimization/model_interface.py b/flowers/optimization/model_interface.py new file mode 100644 index 00000000..6a5970ae --- /dev/null +++ b/flowers/optimization/model_interface.py @@ -0,0 +1,497 @@ +import floris.tools as wfct +import numpy as np +from scipy.interpolate import NearestNDInterpolator +import time + +import flowers.flowers_model as flow +import flowers.optimization.optimization_interface as opt +import flowers.tools.tools as tl + +### THE FOLLOWING CODE IS LEGACY FLOWERS CODE AND IS NOT ACTIVELY MAINTAINED OR +### TESTED IN THE MOST RECENT ARD DEPLOYMENT OF THE FLOWERS METHOD + + +class AEPInterface: + """ + AEPInterface is a high-level user interface to compare AEP estimates between + the FLOWERS (analytical) AEP model and the Conventional (numerical) AEP model. + + Args: + wind_rose (pandas.DataFrame): A dataframe for the wind rose in the FLORIS + format containing the following information: + - 'ws' (float): wind speeds [m/s] + - 'wd' (float): wind directions [deg] + - 'freq_val' (float): frequency for each wind speed and direction + layout_x (numpy.array(float)): x-positions of each turbine [m] + layout_y (numpy.array(float)): y-positions of each turbine [m] + num_terms (int, optional): number of Fourier modes for the FLOWERS model + k (float, optional): wake expansion rate in the FLOWERS model + conventional_model (str, optional): underlying wake model: + - 'jensen' (default) + - 'gauss' + turbine (str, optional): turbine type: + - 'nrel_5MW' (default) + + """ + + ########################################################################### + # Initialization tools + ########################################################################### + + def __init__( + self, + wind_rose, + layout_x, + layout_y, + num_terms=0, + k=0.05, + conventional_model=None, + turbine=None, + ): + + self._wind_rose = wind_rose + self._model = conventional_model + + # Initialize FLOWERS + self.flowers_interface = flow.FlowersInterface( + wind_rose, layout_x, layout_y, num_terms=num_terms, k=k, turbine=turbine + ) + + wd_array = np.array(wind_rose["wd"].unique(), dtype=float) + ws_array = np.array(wind_rose["ws"].unique(), dtype=float) + wd_grid, ws_grid = np.meshgrid(wd_array, ws_array, indexing="ij") + freq_interp = NearestNDInterpolator( + wind_rose[["wd", "ws"]], wind_rose["freq_val"] + ) + freq = freq_interp(wd_grid, ws_grid) + self._freq_2D = freq / np.sum(freq) + + self.floris_interface.reinitialize( + layout_x=layout_x.flatten(), + layout_y=layout_y.flatten(), + wind_directions=wd_array, + wind_speeds=ws_array, + time_series=False, + ) + + def reinitialize( + self, + wind_rose=None, + layout_x=None, + layout_y=None, + num_terms=None, + wd_resolution=0.0, + ws_avg=False, + ): + + # Reinitialize FLOWERS interface + self.flowers_interface.reinitialize( + wind_rose=wind_rose, + layout_x=layout_x, + layout_y=layout_y, + num_terms=num_terms, + ) + + # Reinitialize FLORIS interface + if wind_rose is not None: + self._wind_rose = wind_rose + + wd_array = np.array(wind_rose["wd"].unique(), dtype=float) + ws_array = np.array(wind_rose["ws"].unique(), dtype=float) + wd_grid, ws_grid = np.meshgrid(wd_array, ws_array, indexing="ij") + freq_interp = NearestNDInterpolator( + wind_rose[["wd", "ws"]], wind_rose["freq_val"] + ) + freq = freq_interp(wd_grid, ws_grid) + self._freq_2D = freq / np.sum(freq) + + self.floris_interface.reinitialize( + wind_directions=wd_array, wind_speeds=ws_array, time_series=False + ) + + if layout_x is not None and layout_y is not None: + self.floris_interface.reinitialize( + layout_x=layout_x.flatten(), + layout_y=layout_y.flatten(), + time_series=(np.shape(self._freq_2D)[1] == 1), + ) + elif layout_x is not None and layout_y is None: + self.floris_interface.reinitialize( + layout_x=layout_x.flatten(), + time_series=(np.shape(self._freq_2D)[1] == 1), + ) + elif layout_x is None and layout_y is not None: + self.floris_interface.reinitialize( + layout_y=layout_y.flatten(), + time_series=(np.shape(self._freq_2D)[1] == 1), + ) + + if wd_resolution > 0.0 or ws_avg: + if wd_resolution > 1.0: + wr = tl.resample_wind_direction( + self._wind_rose, wd=np.arange(0, 360, wd_resolution) + ) + else: + wr = self._wind_rose + + if ws_avg: + wr = tl.resample_average_ws_by_wd(wr) + freq = wr.freq_val.to_numpy() + freq /= np.sum(freq) + self._freq_2D = np.expand_dims(freq, 1) + self.floris_interface.reinitialize( + wind_directions=wr.wd, wind_speeds=wr.ws, time_series=True + ) + else: + wr = tl.resample_wind_speed(wr, ws=np.arange(1.0, 26.0, 1.0)) + wd_array = np.array(wr["wd"].unique(), dtype=float) + ws_array = np.array(wr["ws"].unique(), dtype=float) + wd_grid, ws_grid = np.meshgrid(wd_array, ws_array, indexing="ij") + freq_interp = NearestNDInterpolator(wr[["wd", "ws"]], wr["freq_val"]) + freq = freq_interp(wd_grid, ws_grid) + self._freq_2D = freq / np.sum(freq) + self.floris_interface.reinitialize( + wind_directions=wd_array, wind_speeds=ws_array, time_series=False + ) + + ########################################################################### + # AEP methods + ########################################################################### + + def compute_flowers_aep(self, timer=False): + """ + Compute farm AEP using the FLOWERS model. + + Args: + timer (bool, optional): indicate whether wall timer should be + used for calculation. + + Returns: + aep (float): farm AEP [Wh] + elapsed (float, optional): wall time of AEP calculation [s] + + """ + + # Time AEP calculation + if timer: + elapsed = 0 + for _ in range(5): + t = time.time() + aep = self.flowers_interface.calculate_aep() + elapsed += time.time() - t + elapsed /= 5 + return aep, elapsed + else: + aep = self.flowers_interface.calculate_aep() + return aep + + def compute_floris_aep(self, timer=False): + """ + Compute farm AEP using the FLORIS model. + + Args: + timer (bool, optional): indicate whether wall timer should be + used for calculation. + + Returns: + aep (float): farm AEP [Wh] + elapsed (float, optional): wall time of AEP calculation [s] + + """ + + # Time AEP calculation + if timer: + elapsed = 0 + for _ in range(5): + t = time.time() + self.floris_interface.calculate_wake() + aep = np.sum( + self.floris_interface.get_farm_power() * self._freq_2D * 8760 + ) + elapsed += time.time() - t + elapsed /= 5 + return aep, elapsed + + else: + self.floris_interface.calculate_wake() + aep = np.sum(self.floris_interface.get_farm_power() * self._freq_2D * 8760) + return aep + + def compare_aep(self, timer=True, display=True): + """ + Compute farm AEP using both models and compare. The calculation is + repeated an optional number of instances to average computation time. + A table of relevant information is printed to the terminal. + + Args: + iter (int, optional): the number of times AEP should be computed + to average the wall time of each calculation. + num_terms (int, optional): for FLOWERS, the number of Fourier modes + to compute AEP in the range [1, ceiling(num_wind_directions/2)] + ws_avg (bool, optional): for FLORIS, to indicate whether wind speed + should be averaged for each wind direction bin. + wd_resolution (float, optional): for FLORIS, the width of the discrete + wind direction bins to compute AEP + + """ + + if timer: + aep_flowers, time_flowers = self.compute_flowers_aep(timer=True) + aep_floris, time_floris = self.compute_floris_aep(timer=True) + else: + aep_flowers = self.compute_flowers_aep(timer=False) + aep_floris = self.compute_floris_aep(timer=False) + + if display: + print("============================") + print(" AEP Results ") + print(" FLORIS Model: " + str(self._model).capitalize()) + print( + " Number of Turbines: {:.0f}".format( + len(self.flowers_interface.get_layout()[0]) + ) + ) + print( + " FLOWERS Terms: {:.0f}".format( + self.flowers_interface.get_num_modes() + ) + ) + print( + " FLORIS Bins: [{:.0f},{:.0f}]".format( + len(self._freq_2D[:, 0]), len(self._freq_2D[0, :]) + ) + ) + print("----------------------------") + print("FLOWERS AEP: {:.3f} GWh".format(aep_flowers / 1.0e9)) + print("FLORIS AEP: {:.3f} GWh".format(aep_floris / 1.0e9)) + print( + "Percent Difference: {:.1f}%".format( + (aep_flowers - aep_floris) / aep_floris * 100 + ) + ) + if timer: + print("FLOWERS Time: {:.3f} s".format(time_flowers)) + print("FLORIS Time: {:.3f} s".format(time_floris)) + print("Factor of Improvement: {:.1f}x".format(time_floris / time_flowers)) + print("============================") + + if timer: + return (aep_flowers, aep_floris), (time_flowers, time_floris) + else: + return (aep_flowers, aep_floris) + + +class WPLOInterface: + """ + WPLOInterface is a high-level user interface to initialize and run wind plant + layout optimization studies with the FLOWERS (analytical) AEP model and the + Conventional (numerical) AEP model as the objective function. + + Args: + wind_rose (pandas.DataFrame): A dataframe for the wind rose in the FLORIS + format containing the following information: + - 'ws' (float): wind speeds [m/s] + - 'wd' (float): wind directions [deg] + - 'freq_val' (float): frequency for each wind speed and direction + layout_x (numpy.array(float)): x-positions of each turbine [m] + layout_y (numpy.array(float)): y-positions of each turbine [m] + boundaries (list(tuple(float, float))): (x,y) position of each boundary point [m] + num_terms (int, optional): number of Fourier modes for the FLOWERS model + k (float, optional): wake expansion rate in the FLOWERS model + conventional_model (str, optional): underlying wake model: + - 'jensen' (default) + - 'gauss' + turbine (str, optional): turbine type: + - 'nrel_5MW' (default) + + """ + + def __init__( + self, + wind_rose, + layout_x, + layout_y, + boundaries, + num_terms=10, + k=0.05, + conventional_model=None, + turbine=None, + ): + + self._initial_x = layout_x + self._initial_y = layout_y + self._model = conventional_model + self._boundaries = boundaries + + if conventional_model is None or conventional_model == "gauss": + self.floris_interface = wfct.floris_interface.FlorisInterface( + "./input/gauss.yaml" + ) + elif conventional_model == "jensen": + self.floris_interface = wfct.floris_interface.FlorisInterface( + "./input/jensen.yaml" + ) + + # Initialize FLOWERS interface + self.flowers_interface = flow.FlowersInterface( + wind_rose, layout_x, layout_y, num_terms=num_terms, k=k, turbine=turbine + ) + + # Initialize FLORIS interface + wr = tl.resample_wind_direction(wind_rose, wd=np.arange(0, 360, 5.0)) + wr = tl.resample_average_ws_by_wd(wr) + freq = wr.freq_val.to_numpy() + freq /= np.sum(freq) + self._freq_1D = np.expand_dims(freq, 1) + self.floris_interface.reinitialize( + wind_directions=wr.wd, + wind_speeds=wr.ws, + layout_x=layout_x.flatten(), + layout_y=layout_y.flatten(), + time_series=True, + ) + + # Initialize post-processing interface + self.post_processing = wfct.floris_interface.FlorisInterface( + "./input/post.yaml" + ) + wind_rose = tl.resample_wind_speed(wind_rose, ws=np.arange(1.0, 26.0, 1.0)) + wd_array = np.array(wind_rose["wd"].unique(), dtype=float) + ws_array = np.array(wind_rose["ws"].unique(), dtype=float) + wd_grid, ws_grid = np.meshgrid(wd_array, ws_array, indexing="ij") + freq_interp = NearestNDInterpolator( + wind_rose[["wd", "ws"]], wind_rose["freq_val"] + ) + freq = freq_interp(wd_grid, ws_grid) + self._freq_2D = freq / np.sum(freq) + self.post_processing.reinitialize( + layout_x=layout_x.flatten(), + layout_y=layout_y.flatten(), + wind_directions=wd_array, + wind_speeds=ws_array, + time_series=False, + ) + + # Calculate initial AEP + self.post_processing.calculate_wake() + self._aep_initial = np.sum( + self.post_processing.get_farm_power() * self._freq_2D * 8760 + ) + + def run_optimization( + self, + optimizer, + gradient="analytic", + solver="SLSQP", + scale=1e3, + tol=1e-2, + timer=None, + history="hist.hist", + output="out.out", + ): + """ + Run a Wind Plant Layout Optimization study with either the FLOWERS + or Conventional optimizer. + + Args: + optimizer (str): the objective function to use in the study: + - "flowers" + - "conventional" + solver (str, optional): the optimization algorithm to use: + - "SLSQP" (default) + - "SNOPT" + timer (int, optional): time limit [s] + history (str, optional): file name for pyoptsparse history file + output (str, optional): file name for solver output file + + Returns: + solution (dict): relevant information from the optimization solution: + - "init_x" (numpy.array(float)): initial x-positions of each turbine [m] + - "init_y" (numpy.array(float)): initial y-positions of each turbine [m] + - "init_aep" (float): initial plant AEP [Wh] + - "opt_x" (numpy.array(float)): optimized x-positions of each turbine [m] + - "opt_y" (numpy.array(float)): optimized y-positions of each turbine [m] + - "opt_aep" (float): optimized plant AEP [Wh] + - "hist_x" (numpy.array(float,float)): x-positions of each turbine at each solver iteration [m] + - "hist_y" (numpy.array(float,float)): y-positions of each turbine at each solver iteration [m] + - "hist_aep" (numpy.array(float)): plant AEP at each solver iteration [Wh] + - "iter" (int): number of major iterations taken by the solver + - "obj_calls" (int): number of AEP function evaluations + - "grad_calls" (int): number of gradient evaluations + - "total_time" (float): total solve time + - "obj_time" (float): time spent evaluating objective function + - "grad_time" (float): time spent evaluating gradients + - "solver_time" (float): time spent solving optimization problem + + """ + + # Instantiate optimizer class with user inputs + if optimizer == "flowers": + prob = opt.FlowersOptimizer( + self.flowers_interface, + self._initial_x, + self._initial_y, + self._boundaries, + grad=gradient, + solver=solver, + scale=scale, + tol=tol, + timer=timer, + history_file=history, + output_file=output, + ) + elif optimizer == "conventional": + prob = opt.ConventionalOptimizer( + self.floris_interface, + self._freq_1D, + self._initial_x, + self._initial_y, + self._boundaries, + grad=gradient, + solver=solver, + scale=scale, + tol=tol, + timer=timer, + history_file=history, + output_file=output, + ) + + # Solve optimization problem + print("Solving layout optimization problem.") + sol = prob.optimize() + print("Optimization complete: " + str(sol.optInform["text"])) + + # Define solution dictionary and gather data + self.solution = dict() + self.solution["init_x"] = self._initial_x + self.solution["init_y"] = self._initial_y + self.solution["opt_x"], self.solution["opt_y"] = prob.parse_sol_vars(sol) + self.solution["total_time"] = float(sol.optTime) + self.solution["obj_time"] = float(sol.userObjTime) + self.solution["grad_time"] = float(sol.userSensTime) + self.solution["solver_time"] = float(sol.optCodeTime) + self.solution["obj_calls"] = int(sol.userObjCalls) + self.solution["grad_calls"] = int(sol.userSensCalls) + self.solution["exit_code"] = sol.optInform["text"] + self.solution["init_aep"] = self._aep_initial + + # Get number of iterations + with open(output, "r") as fp: + for line in fp: + if "No. of major iterations" in line: + self.solution["iter"] = int(line.split()[4]) + break + + # Compute optimal AEP + self.post_processing.reinitialize( + layout_x=self.solution["opt_x"].flatten(), + layout_y=self.solution["opt_y"].flatten(), + time_series=False, + ) + self.post_processing.calculate_wake() + self._aep_final = np.sum( + self.post_processing.get_farm_power() * self._freq_2D * 8760 + ) + self.solution["opt_aep"] = self._aep_final + + return self.solution diff --git a/flowers/optimization/optimization_interface.py b/flowers/optimization/optimization_interface.py new file mode 100644 index 00000000..8e30d571 --- /dev/null +++ b/flowers/optimization/optimization_interface.py @@ -0,0 +1,404 @@ +import numpy as np +import pyoptsparse + +### THE FOLLOWING CODE IS LEGACY FLOWERS CODE AND IS NOT ACTIVELY MAINTAINED OR +### TESTED IN THE MOST RECENT ARD DEPLOYMENT OF THE FLOWERS METHOD + + +class LayoutOptimizer: + """ + LayoutOptimizer is a base class for wind plant layout optimization, + acting as a wrapper for pyOptSparse. + + """ + + def _base_init_( + self, + layout_x, + layout_y, + boundaries, + solver="SNOPT", + tol=1e-2, + timer=None, + options=None, + history_file="hist.hist", + output_file="out.out", + ): + + # Save boundary information + self._boundaries = np.array(boundaries).T + self._nbounds = len(self._boundaries[0]) + + # Compute edge information + self._boundary_edge = np.roll(self._boundaries, -1, axis=1) - self._boundaries + self._boundary_len = np.sqrt( + self._boundary_edge[0] ** 2 + self._boundary_edge[1] ** 2 + ) + self._boundary_norm = ( + np.array([self._boundary_edge[1], -self._boundary_edge[0]]) + / self._boundary_len + ) + self._boundary_int = ( + np.roll(self._boundary_norm, 1, axis=1) + self._boundary_norm + ) / 2 + + # Position normalization + self._xmin = np.min(self._boundaries[0]) + self._xmax = np.max(self._boundaries[0]) + self._ymin = np.min(self._boundaries[1]) + self._ymax = np.max(self._boundaries[1]) + + self._x0 = layout_x + self._y0 = layout_y + self._nturbs = len(layout_x) + + # Optimization initialization + self.solver = solver + # self.storeHistory = history_file + self.timeLimit = timer + self.optProb = pyoptsparse.Optimization("layout", self._obj_func) + + self.optProb = self.add_var_group(self.optProb) + self.optProb = self.add_con_group(self.optProb) + self.optProb.addObj("obj") + + # Optimizer options + if options is not None: + self.optOptions = options + elif solver == "SNOPT": + self.optOptions = { + # "Print file": output_file, + # "iSumm": 0, + "Summary file": output_file, + "iPrint": 0, + "Major optimality tolerance": tol, + "Minor optimality tolerance": 1e-4, + "Major feasibility tolerance": 1e-4, + "Minor feasibility tolerance": 1e-4, + "Scale option": 0, + # "Verify level": 3, + } + elif solver == "SLSQP": + self.optOptions = { + "ACC": 1e-6, + "IFILE": output_file, + "MAXIT": 100, + } + elif solver == "NSGA2": + self.optOptions = { + "maxGen": 10, + } + + exec("self.opt = pyoptsparse." + self.solver + "(options=self.optOptions)") + + def _norm(self, val, x1, x2): + """Method to normalize turbine positions""" + return (val - x1) / (x2 - x1) + + def _unnorm(self, val, x1, x2): + """Method to dimensionalize turbine positions""" + return np.array(val) * (x2 - x1) + x1 + + def optimize(self): + """Method to initiate optimization.""" + self._optimize() + return self.sol + + ########################################################################### + # User constraint function + ########################################################################### + def _boundary_constraint(self, gradient=False): + # Transform inputs + points = np.array([self._x, self._y]) + + # Compute distances from turbines to boundary points + a = np.zeros((self._nturbs, 2, self._nbounds)) + for i in range(self._nturbs): + a[i] = np.expand_dims(points[:, i].T, axis=-1) - self._boundaries + + # Compute projections + a_edge = np.sum(a * self._boundary_edge, axis=1) / self._boundary_len + a_int = np.sum(a * self._boundary_norm, axis=1) + sigma = np.sign(np.sum(a * self._boundary_int, axis=1)) + + # Initialize signed distance containers + C = np.zeros(self._nturbs) + D = np.zeros(self._nbounds) + if gradient: + Cx = np.zeros(self._nturbs) + Cy = np.zeros(self._nturbs) + + # Compute signed distance + for i in range(self._nturbs): + for k in range(self._nbounds): + if a_edge[i, k] < 0: + D[k] = np.sqrt(a[i, 0, k] ** 2 + a[i, 1, k] ** 2) * sigma[i, k] + elif a_edge[i, k] > self._boundary_len[k]: + D[k] = ( + np.sqrt( + a[i, 0, (k + 1) % self._nbounds] ** 2 + + a[i, 1, (k + 1) % self._nbounds] ** 2 + ) + * sigma[i, (k + 1) % self._nbounds] + ) + else: + D[k] = a_int[i, k] + + # Select minimum distance + idx = np.argmin(np.abs(D)) + C[i] = D[idx] + + if gradient: + if a_edge[i, idx] < 0: + Cx[i] = (points[0, i] - self._boundaries[0, idx]) / np.sqrt( + (self._boundaries[0, idx] - points[0, i]) ** 2 + + (self._boundaries[1, idx] - points[1, i]) ** 2 + ) + Cy[i] = (points[1, i] - self._boundaries[1, idx]) / np.sqrt( + (self._boundaries[0, idx] - points[0, i]) ** 2 + + (self._boundaries[1, idx] - points[1, i]) ** 2 + ) + elif a_edge[i, idx] > self._boundary_len[idx]: + Cx[i] = ( + points[0, i] - self._boundaries[0, (idx + 1) % self._nbounds] + ) / np.sqrt( + (self._boundaries[0, (idx + 1) % self._nbounds] - points[0, i]) + ** 2 + + ( + self._boundaries[1, (idx + 1) % self._nbounds] + - points[1, i] + ) + ** 2 + ) + Cy[i] = ( + points[1, i] - self._boundaries[1, (idx + 1) % self._nbounds] + ) / np.sqrt( + (self._boundaries[0, (idx + 1) % self._nbounds] - points[0, i]) + ** 2 + + ( + self._boundaries[1, (idx + 1) % self._nbounds] + - points[1, i] + ) + ** 2 + ) + else: + Cx[i] = ( + self._boundaries[1, (idx + 1) % self._nbounds] + - self._boundaries[1, idx] + ) / self._boundary_len[idx] + Cy[i] = ( + self._boundaries[0, idx] + - self._boundaries[0, (idx + 1) % self._nbounds] + ) / self._boundary_len[idx] + + if gradient: + return C, Cx, Cy + else: + return C + + ########################################################################### + # pyOptSparse wrapper functions + ########################################################################### + def _optimize(self): + if self.gradient: + if self.timeLimit is not None: + self.sol = self.opt( + self.optProb, sens=self._sens_func + ) # , timeLimit=self.timeLimit) #storeHistory=self.storeHistory + else: + self.sol = self.opt(self.optProb, sens=self._sens_func) + else: + if self.timeLimit is not None: + self.sol = self.opt(self.optProb, timeLimit=self.timeLimit) + else: + self.sol = self.opt(self.optProb, sens="FD") + + def parse_opt_vars(self, varDict): + # self._x = self._unnorm(varDict["x"], self._xmin, self._xmax) + # self._y = self._unnorm(varDict["y"], self._ymin, self._ymax) + self._x = varDict["x"] + self._y = varDict["y"] + + def parse_sol_vars(self, sol): + # return np.array(self._unnorm(sol.getDVs()["x"], self._xmin, self._xmax)), np.array(self._unnorm(sol.getDVs()["y"], self._ymin, self._ymax)) + return np.array(sol.getDVs()["x"]), np.array(sol.getDVs()["y"]) + + def parse_hist_vars(self, hist): + val = hist.getValues(names=["x", "y"], major=True) + # return np.array(self._unnorm(val['x'], self._xmin, self._xmax)), np.array(self._unnorm(val['y'], self._ymin, self._ymax)) + return np.array(val["x"]), np.array(val["y"]) + + def add_var_group(self, optProb): + optProb.addVarGroup("x", self._nturbs, type="c", value=self._x0) + optProb.addVarGroup("y", self._nturbs, type="c", value=self._y0) + return optProb + + def add_con_group(self, optProb): + optProb.addConGroup("con", self._nturbs, upper=0.0) + return optProb + + def compute_cons(self, funcs): + funcs["con"] = self._boundary_constraint() + return funcs + + +class FlowersOptimizer(LayoutOptimizer): + """ + Child class of LayoutOptimizer for the FLOWERS-based layout optimizer. + + Args: + flowers_interface (FlowersInterface): FLOWERS interface for calculating AEP + layout_x (numpy.array(float)): x-positions of each turbine [m] + layout_y (numpy.array(float)): y-positions of each turbine [m] + boundaries (list(tuple(float, float))): (x,y) position of each boundary point [m] + solver (str, optional): the optimization algorithm to use: + - "SLSQP" (default) + - "SNOPT" + timer (int, optional): time limit [s] + history (str, optional): file name for pyoptsparse history file + output (str, optional): file name for solver output file + + """ + + def __init__( + self, + flowers_interface, + layout_x, + layout_y, + boundaries, + grad="analytical", + solver="SNOPT", + scale=1e3, + tol=1e-2, + timer=None, + history_file="hist.hist", + output_file="out.out", + ): + self.model = flowers_interface + if grad == "analytical": + self.gradient = True + elif grad == "numerical": + self.gradient = False + aep_initial = self.model.calculate_aep(gradient=False) + self._scale = aep_initial / scale + self._base_init_( + layout_x, + layout_y, + boundaries, + solver=solver, + tol=tol, + timer=timer, + history_file=history_file, + output_file=output_file, + ) + + def _obj_func(self, varDict): + # Parse the variable dictionary + self.parse_opt_vars(varDict) + + # Update turbine map with turbince locations + self.model.reinitialize(layout_x=self._x, layout_y=self._y) + + # Compute the objective function + funcs = {} + funcs["obj"] = -1 * self.model.calculate_aep() / self._scale + + # Compute constraints, if any are defined for the optimization + funcs = self.compute_cons(funcs) + + fail = False + return funcs, fail + + # Optionally, the user can supply the optimization with gradients + def _sens_func(self, varDict, funcs): + # Parse the variable dictionary + self.parse_opt_vars(varDict) + + self.model.reinitialize(layout_x=self._x, layout_y=self._y) + + _, tmp = self.model.calculate_aep(gradient=True) + funcsSens = {} + funcsSens["obj"] = { + "x": -tmp[:, 0] / self._scale, + "y": -tmp[:, 1] / self._scale, + } + + _, tmpx, tmpy = self._boundary_constraint(gradient=True) + funcsSens["con"] = {"x": np.diag(tmpx), "y": np.diag(tmpy)} + + fail = False + + return funcsSens, fail + + +class ConventionalOptimizer(LayoutOptimizer): + """ + Child class of LayoutOptimizer for the FLORIS-based layout optimizer. + + Args: + floris_interface (FlorisInterface): FLORIS interface for calculating AEP + layout_x (numpy.array(float)): x-positions of each turbine [m] + layout_y (numpy.array(float)): y-positions of each turbine [m] + boundaries (list(tuple(float, float))): (x,y) position of each boundary point [m] + solver (str, optional): the optimization algorithm to use: + - "SLSQP" (default) + - "SNOPT" + timer (int, optional): time limit [s] + history (str, optional): file name for pyoptsparse history file + output (str, optional): file name for solver output file + + """ + + def __init__( + self, + floris_interface, + freq_val, + layout_x, + layout_y, + boundaries, + grad="analytical", + solver="SNOPT", + scale=1e3, + tol=1e-2, + timer=None, + history_file="hist.hist", + output_file="out.out", + ): + self.model = floris_interface + self.gradient = False + self._freq_1D = freq_val + self.model.calculate_wake() + self._scale = np.sum(self.model.get_farm_power() * self._freq_1D * 8760) / scale + self._base_init_( + layout_x, + layout_y, + boundaries, + solver=solver, + timer=timer, + history_file=history_file, + output_file=output_file, + ) + + def _obj_func(self, varDict): + # Parse the variable dictionary + self.parse_opt_vars(varDict) + + # Update turbine map with turbince locations + self.model.reinitialize( + layout_x=self._x.flatten(), layout_y=self._y.flatten(), time_series=True + ) + + # Compute the objective function + self.model.calculate_wake() + funcs = {} + funcs["obj"] = ( + -1 + * np.sum(self.model.get_farm_power() * self._freq_1D * 8760) + / self._scale + ) + + # Compute constraints, if any are defined for the optimization + funcs = self.compute_cons(funcs) + + fail = False + return funcs, fail diff --git a/flowers/tools.py b/flowers/tools.py new file mode 100644 index 00000000..443e569b --- /dev/null +++ b/flowers/tools.py @@ -0,0 +1,652 @@ +import copy + +import numpy as np +import pandas as pd +import pickle +from shapely.geometry import Polygon, Point + +# import floris.tools.wind_rose as rose + + +########################################################################### +# Layout generation +########################################################################### + + +def random_layout(boundaries=[], n_turb=0, D=126.0, min_dist=2.0, seed_val=None): + """ + Generate a random wind farm layout within the specified boundaries. + Minimum spacing between turbines is 2D. + + Args: + boundaries (list(tuple)): boundary vertices in the form + [(x0,y0), (x1,y1), ... , (xN,yN)] + n_turb (int): number of turbines + D (float): rotor diameter [m] + min_dist (float): enforced minimum spacing between turbine centers + normalized by rotor diameter + seed_val (int, optional): random number generator seed + + Args: + xx (np.array): x-positions of each turbine + yy (np.array): y-positions of each turbine + + """ + + print("Generating wind farm layout.") + + # Verify that boundaries and turbines are supplied + if not boundaries: + raise ValueError("Must supply boundaries to generate wind farm.") + + if n_turb <= 0: + raise ValueError("Must supply number of turbines.") + + # Initialize RNG and containers + if seed_val != None: + np.random.seed(seed_val) + + xx = -1000 * np.ones(n_turb) + yy = -1000 * np.ones(n_turb) + + xmin = np.min([tup[0] for tup in boundaries]) + xmax = np.max([tup[0] for tup in boundaries]) + ymin = np.min([tup[1] for tup in boundaries]) + ymax = np.max([tup[1] for tup in boundaries]) + + # Generate boundary polygon + poly = Polygon(boundaries) + + # Generate new positions and check minimum spacing and boundary + for i in range(n_turb): + prop_x = np.random.uniform(low=xmin, high=xmax) + prop_y = np.random.uniform(low=ymin, high=ymax) + pt = Point(prop_x, prop_y) + while np.any( + np.sqrt((prop_x - xx) ** 2 + (prop_y - yy) ** 2) < min_dist * D + ) or not pt.within(poly): + prop_x = np.random.uniform(low=xmin, high=xmax) + prop_y = np.random.uniform(low=ymin, high=ymax) + pt = Point(prop_x, prop_y) + xx[i] = prop_x + yy[i] = prop_y + + return xx, yy + + +def discrete_layout(n_turb=0, D=126.0, min_dist=3.0, seed_val=None, spacing=False): + """ + Generate a random wind farm layout within the specified boundaries. + Minimum spacing between turbines is 2D. + + Args: + boundaries (list(tuple)): boundary vertices in the form + [(x0,y0), (x1,y1), ... , (xN,yN)] + n_turb (int): number of turbines + D (float): rotor diameter [m] + min_dist (float): enforced minimum spacing between turbine centers + normalized by rotor diameter + idx (int, optional): random number generator seed + + Args: + xx (np.array): x-positions of each turbine + yy (np.array): y-positions of each turbine + + """ + + print("Generating wind farm layout.") + + if n_turb <= 0: + raise ValueError("Must supply number of turbines.") + + # Initialize RNG and containers + if seed_val != None: + np.random.seed(seed_val) + + # Indices of discrete grid + sx = n_turb + 1 + sy = 6 + x_idx = np.random.randint(0, sx, n_turb) + y_idx = np.random.randint(0, sy, n_turb) + pts = [(x_idx[i], y_idx[i]) for i in range(n_turb)] + while len(np.unique(pts)) < len(pts): + tmp = np.unique(pts) + new_set = [] + for i in range(n_turb): + if i not in tmp: + new_set.append(i) + x_idx[new_set] = np.random.randint(0, sx, len(new_set)) + y_idx[new_set] = np.random.randint(0, sy, len(new_set)) + pts = [(x_idx[i], y_idx[i]) for i in range(n_turb)] + + # Check that all combinations of x,y are unique + + xx = np.array(min_dist * D * x_idx) + yy = np.array(min_dist * D * y_idx) + + if spacing: + x_rel = (xx - np.reshape(xx, (-1, 1))) / D + y_rel = (yy - np.reshape(yy, (-1, 1))) / D + r_rel = np.sqrt(x_rel**2 + y_rel**2) + r_rel = np.ma.masked_where(np.eye(len(xx)), r_rel) + ss = np.mean(np.min(r_rel, -1)) + return xx, yy, ss + else: + return xx, yy + + +def load_layout(idx, case, boundaries=True): + file = "./layouts/" + case + str(idx) + ".p" + with open(file, "rb") as infile: + layout_x, layout_y, boundaries_val = pickle.load(infile) + + if boundaries: + return layout_x, layout_y, boundaries_val + else: + return layout_x, layout_y + + +########################################################################### +# Wind rose sampling +########################################################################### + + +def load_wind_rose(idx): + """ + Load a locally-stored wind rose saved to a pickle file. + See show_wind_roses.py to visualize all wind rose options. + + Args: + idx (int): index of desired wind rose + + Returns: + df (pandas.DataFrame): A dataframe for the wind rose in the FLORIS + format containing the following information: + - 'ws' (float): wind speeds [m/s] + - 'wd' (float): wind directions [deg] + - 'freq_val' (float): frequency for each wind speed and direction + + """ + + print("Generating wind rose.") + file_name = "./wind_roses/wr" + str(idx) + ".p" + df = pd.read_pickle(file_name) + + return df + + +def resample_wind_direction(df, wd=np.arange(0, 360, 5.0)): + """ + Resample wind direction bins using new specified bin center values. + (Copied from FLORIS) + + Args: + df (pandas.DataFrame): Wind rose DataFrame containing the following + columns: + - 'wd': Wind direction bin center values (deg). + - 'ws': Wind speed bin center values (m/s). + - 'freq_val': The frequency of occurence of the + wind conditions in the other columns. + + wd (np.array, optional): List of new wind direction center bins + (deg). Defaults to np.arange(0, 360, 5.). + + Returns: + New wind rose DataFrame containing the following columns: + - 'wd': New wind direction bin center values from wd argument (deg). + - 'ws': Resampled wind speed bin center values (m/s). + - 'freq_val': The resampled frequency of occurence of the + wind conditions in the other columns. + """ + + # Make a copy of incoming dataframe + df = df.copy(deep=True) + + # Get the wind step + wd_step = wd[1] - wd[0] + + # Get bin edges + wd_edges = wd - wd_step / 2.0 + wd_edges = np.append(wd_edges, np.array(wd[-1] + wd_step / 2.0)) + + # Get the overhangs + negative_overhang = wd_edges[0] + positive_overhang = wd_edges[-1] - 360.0 + + # Potentially wrap high angle direction to negative for correct binning + tmp = df.wd + tmp = np.where(tmp < 0.0, tmp + 360.0, tmp) + tmp = np.where(tmp >= 360.0, tmp - 360.0, tmp) + df["wd"] = tmp + + if negative_overhang < 0: + # print("Correcting negative Overhang:%.1f" % negative_overhang) + df["wd"] = np.where( + df.wd.values >= 360.0 + negative_overhang, + df.wd.values - 360.0, + df.wd.values, + ) + + # Check on other side + if positive_overhang > 0: + # print("Correcting positive Overhang:%.1f" % positive_overhang) + df["wd"] = np.where( + df.wd.values <= positive_overhang, df.wd.values + 360.0, df.wd.values + ) + + # Cut into bins + df["wd"] = pd.cut(df.wd, wd_edges, labels=wd) + + # Regroup + df = df.groupby([c for c in df.columns if c != "freq_val"]).sum() + + # Fill nans + df = df.fillna(0) + + # Reset the index + df = df.reset_index() + + # Set to float Re-wrap + for c in [c for c in df.columns if c != "freq_val"]: + df[c] = df[c].astype(float) + df[c] = df[c].astype(float) + + tmp = df.wd + tmp = np.where(tmp < 0.0, tmp + 360.0, tmp) + tmp = np.where(tmp >= 360.0, tmp - 360.0, tmp) + df["wd"] = tmp + + return df + + +def resample_wind_speed(df, ws=np.arange(0, 26, 1.0)): + """ + Resample wind speed bins using new specified bin center values. + (Copied from FLORIS) + + Args: + df (pandas.DataFrame): Wind rose DataFrame containing the following + columns: + - 'wd': Wind direction bin center values (deg). + - 'ws': Wind speed bin center values (m/s). + - 'freq_val': The frequency of occurence of the + wind conditions in the other columns. + + ws (np.array, optional): List of new wind direction center bins + (m/s). Defaults to np.arange(0, 26, 1.0). + + Returns: + New wind rose DataFrame containing the following columns: + - 'wd': New wind direction bin center values from wd argument (deg). + - 'ws': Resampled wind speed bin center values (m/s). + - 'freq_val': The resampled frequency of occurence of the + wind conditions in the other columns. + """ + # Make a copy of incoming dataframe + df = df.copy(deep=True) + + # Get the wind step + ws_step = ws[1] - ws[0] + + # Ws + ws_edges = ws - ws_step / 2.0 + ws_edges = np.append(ws_edges, np.array(ws[-1] + ws_step / 2.0)) + + # Cut wind speed onto bins + df["ws"] = pd.cut(df.ws, ws_edges, labels=ws) + + # Regroup + df = df.groupby([c for c in df.columns if c != "freq_val"]).sum() + + # Fill nans + df = df.fillna(0) + + # Reset the index + df = df.reset_index() + + # Set to float + for c in [c for c in df.columns if c != "freq_val"]: + df[c] = df[c].astype(float) + df[c] = df[c].astype(float) + + return df + + +def resample_average_ws_by_wd(wind_rose): + """ + Calculate the mean wind speed for each wind direction bin + and resample the wind rose. (Copied from FLORIS) + + Args: + wind_rose (WindRose): Wind rose object containing the following + columns: + - 'wd': Wind direction bin center values (deg). + - 'ws': Wind speed bin center values (m/s). + - 'freq_val': The frequency of occurence of the + wind conditions in the other columns. + + Returns: + New wind rose DataFrame containing the following columns: + - 'wd': Wind direction bin center values (deg). + - 'ws': Resampled average wind speed values (m/s). + - 'freq_val': The resampled frequency of occurence of the + wind conditions in the other columns. + + """ + # Make a copy of incoming FLORIS WindRose object + df = copy.deepcopy(wind_rose) + # wind_rose = copy.deepcopy(wind_rose) + + ws_avg = [] + + for val in df.wd.unique(): + ws_avg.append( + np.array( + df.loc[df["wd"] == val]["ws"] * df.loc[df["wd"] == val]["freq_val"] + ).sum() + / df.loc[df["wd"] == val]["freq_val"].sum() + ) + + # freq_avg = np.sum(wind_rose.freq_table, axis=1) + # ws_avg = np.sum( + # wind_rose.freq_table * wind_rose.wind_speeds.reshape(1, -1), axis=1 + # ) / freq_avg + + # print(wind_rose.wind_directions) + # print(wind_rose.freq_table) + # print(ws_avg) + # lkj + # df = pd.DataFrame() + # df["wd"] = wind_rose.wind_directions + # df["ws"] = ws_avg + # df["freq_val"] = freq_avg + + # Regroup + df = df.groupby("wd").sum() + + df["ws"] = ws_avg + + # Reset the index + df = df.reset_index() + + # Set to float + df["ws"] = df.ws.astype(float) + df["wd"] = df.wd.astype(float) + + return df + + +########################################################################### +# Turbine parameter tables +########################################################################### + + +def ct_lookup(u, turbine_type, ct=None): + """ + Look-up table for thrust coefficient of the NREL 5 MW turbine. + + Args: + u (float): normalized inflow wind speed + + Returns: + ct (float): thrust coefficient + + """ + + if ct != None: + ct_table = np.array([0.0, 0.0, ct, ct, 0.0, 0.0]) + u_table = 1 / 25.0 * np.array([0.0, 2.0, 2.5, 25.01, 25.02, 50.0]) + elif turbine_type == "nrel_5MW": + ct_table = np.array( + [ + 0.0, + 0.0, + 0.0, + 0.99, + 0.99, + 0.97373036, + 0.92826162, + 0.89210543, + 0.86100905, + 0.835423, + 0.81237673, + 0.79225789, + 0.77584769, + 0.7629228, + 0.76156073, + 0.76261984, + 0.76169723, + 0.75232027, + 0.74026851, + 0.72987175, + 0.70701647, + 0.54054532, + 0.45509459, + 0.39343381, + 0.34250785, + 0.30487242, + 0.27164979, + 0.24361964, + 0.21973831, + 0.19918151, + 0.18131868, + 0.16537679, + 0.15103727, + 0.13998636, + 0.1289037, + 0.11970413, + 0.11087113, + 0.10339901, + 0.09617888, + 0.09009926, + 0.08395078, + 0.0791188, + 0.07448356, + 0.07050731, + 0.06684119, + 0.06345518, + 0.06032267, + 0.05741999, + 0.05472609, + 0.0, + 0.0, + ] + ) + u_table = ( + 1 + / 25.0 + * np.array( + [ + 0.0, + 2.0, + 2.5, + 3.0, + 3.5, + 4.0, + 4.5, + 5.0, + 5.5, + 6.0, + 6.5, + 7.0, + 7.5, + 8.0, + 8.5, + 9.0, + 9.5, + 10.0, + 10.5, + 11.0, + 11.5, + 12.0, + 12.5, + 13.0, + 13.5, + 14.0, + 14.5, + 15.0, + 15.5, + 16.0, + 16.5, + 17.0, + 17.5, + 18.0, + 18.5, + 19.0, + 19.5, + 20.0, + 20.5, + 21.0, + 21.5, + 22.0, + 22.5, + 23.0, + 23.5, + 24.0, + 24.5, + 25.0, + 25.01, + 25.02, + 50.0, + ] + ) + ) + else: + raise NotImplementedError("your specified turbine type was not found.") + + return np.interp(u, u_table, ct_table) + + +def cp_lookup(u, turbine_type, cp=None): + """ + Look-up table for power coefficient of the NREL 5 MW turbine. + + Args: + u (float): normalized inflow wind speed + + Returns: + cp (float): power coefficient + + """ + if cp != None: + cp_table = np.array([0.0, 0.0, cp, cp, 0.0, 0.0]) + u_table = 1 / 25.0 * np.array([0.0, 2.0, 2.5, 25.01, 25.02, 50.0]) + elif turbine_type == "nrel_5MW": + cp_table = np.array( + [ + 0.0, + 0.0, + 0.0, + 0.178085, + 0.289075, + 0.349022, + 0.384728, + 0.406059, + 0.420228, + 0.428823, + 0.433873, + 0.436223, + 0.436845, + 0.436575, + 0.436511, + 0.436561, + 0.436517, + 0.435903, + 0.434673, + 0.433230, + 0.430466, + 0.378869, + 0.335199, + 0.297991, + 0.266092, + 0.238588, + 0.214748, + 0.193981, + 0.175808, + 0.159835, + 0.145741, + 0.133256, + 0.122157, + 0.112257, + 0.103399, + 0.095449, + 0.088294, + 0.081836, + 0.075993, + 0.070692, + 0.065875, + 0.061484, + 0.057476, + 0.053809, + 0.050447, + 0.047358, + 0.044518, + 0.041900, + 0.039483, + 0.0, + 0.0, + ] + ) + u_table = ( + 1 + / 25.0 + * np.array( + [ + 0.0, + 2.0, + 2.5, + 3.0, + 3.5, + 4.0, + 4.5, + 5.0, + 5.5, + 6.0, + 6.5, + 7.0, + 7.5, + 8.0, + 8.5, + 9.0, + 9.5, + 10.0, + 10.5, + 11.0, + 11.5, + 12.0, + 12.5, + 13.0, + 13.5, + 14.0, + 14.5, + 15.0, + 15.5, + 16.0, + 16.5, + 17.0, + 17.5, + 18.0, + 18.5, + 19.0, + 19.5, + 20.0, + 20.5, + 21.0, + 21.5, + 22.0, + 22.5, + 23.0, + 23.5, + 24.0, + 24.5, + 25.0, + 25.01, + 25.02, + 50.0, + ] + ) + ) + else: + raise NotImplementedError("your specified turbine type was not found.") + + return np.interp(u, u_table, cp_table) diff --git a/flowers/visualization.py b/flowers/visualization.py new file mode 100644 index 00000000..cb6bbb73 --- /dev/null +++ b/flowers/visualization.py @@ -0,0 +1,574 @@ +# FLOWERS + +# Michael LoCascio + +import matplotlib.animation as animation +import matplotlib.cm as cm +import matplotlib.pyplot as plt +import numpy as np + +import flowers.tools as tl + +########################################################################### +# Wind rose methods +########################################################################### + + +def plot_wind_rose( + wind_rose, + ax=None, + color_map="viridis_r", + ws_right_edges=np.array([5, 10, 15, 20, 25]), + wd_bins=np.arange(0, 360, 15.0), + legend_kwargs={}, +): + """ + Plots a wind rose showing the frequency of occurence + of the specified wind direction and wind speed bins. If no axis is + provided, a new one is created. (Copied from FLORIS) + + Args: + ax (:py:class:`matplotlib.pyplot.axes`, optional): The figure axes + on which the wind rose is plotted. Defaults to None. + color_map (str, optional): Colormap to use. Defaults to 'viridis_r'. + ws_right_edges (np.array, optional): The upper bounds of the wind + speed bins (m/s). The first bin begins at 0. Defaults to + np.array([5, 10, 15, 20, 25]). + wd_bins (np.array, optional): The wind direction bin centers used + for plotting (deg). Defaults to np.arange(0, 360, 15.). + legend_kwargs (dict, optional): Keyword arguments to be passed to + ax.legend(). + + Returns: + :py:class:`matplotlib.pyplot.axes`: A figure axes object containing + the plotted wind rose. + + """ + # Resample data onto bins + wr = wind_rose.copy(deep=True) + df_plot = tl.resample_wind_direction(wr, wd=wd_bins) + + # Make labels for wind speed based on edges + ws_step = ws_right_edges[1] - ws_right_edges[0] + ws_labels = ["%d-%d m/s" % (w - ws_step, w) for w in ws_right_edges] + + # Grab the wd_step + wd_step = wd_bins[1] - wd_bins[0] + + # Set up figure + if ax is None: + _, ax = plt.subplots(subplot_kw=dict(polar=True)) + + # Get a color array + color_array = cm.get_cmap(color_map, len(ws_right_edges)) + + for wd_idx, wd in enumerate(wd_bins): + rects = list() + df_plot_sub = df_plot[df_plot.wd == wd] + for ws_idx, ws in enumerate(ws_right_edges[::-1]): + plot_val = df_plot_sub[ + df_plot_sub.ws <= ws + ].freq_val.sum() # Get the sum of frequency up to this wind speed + rects.append( + ax.bar( + np.radians(wd), + plot_val, + width=0.9 * np.radians(wd_step), + color=color_array(ws_idx), + edgecolor="k", + linewidth=0.5, + label=ws_labels[ws_idx], + ) + ) + # break + + # Configure the plot + ax.legend( + reversed(rects), + ws_labels, + loc="lower left", + bbox_to_anchor=(0.55 + np.cos(0.55) / 2, 0.4 + np.sin(0.55) / 2), + **legend_kwargs, + ) + ax.set_theta_direction(-1) + ax.set_theta_offset(np.pi / 2.0) + ax.set_theta_zero_location("N") + ax.tick_params(axis="x", which="major", pad=-1) # 1 for AEP paper + ax.set_xticklabels(["N", "NE", "E", "SE", "S", "SW", "W", "NW"]) + ax.set_yticklabels([]) + ax.set_axisbelow(True) + + return ax + + +########################################################################### +# Layout methods +########################################################################### + + +def plot_layout( + layout_x, layout_y, D=126.0, boundaries=None, norm=True, ax=None, color="tab:blue" +): + """ + Plot a wind farm layout. The turbine markers are properly scaled + relative to the domain. + + Args: + layout_x (numpy.array(float)): x-positions of each turbine [m] + layout_y (numpy.array(float)): y-positions of each turbine [m] + D (float): rotor diameter [m] + norm (bool): dictates whether the plot should be scaled by + rotor diameter. Defaults to True. + ax (:py:class:`matplotlib.pyplot.axes`, optional): axis to + plot layout + + Returns: + ax (:py:class:`matplotlib.pyplot.axes`, optional): axis + after plotting and formatting + + """ + + if ax is None: + _, ax = plt.subplots() + + if norm: + xx = layout_x / D + yy = layout_y / D + r = 0.5 + xlab = "x/D" + ylab = "y/D" + + else: + xx = layout_x + yy = layout_x + r = D / 2 + xlab = "x [m]" + ylab = "y [m]" + + if boundaries is not None: + verts = np.array(boundaries) / D + for i in range(len(verts)): + if i == len(verts) - 1: + ax.plot([verts[i][0], verts[0][0]], [verts[i][1], verts[0][1]], "black") + else: + ax.plot( + [verts[i][0], verts[i + 1][0]], + [verts[i][1], verts[i + 1][1]], + "black", + ) + + ax.scatter(xx, yy, s=0.01) + for x, y in zip(xx, yy): + ax.add_patch(plt.Circle((x, y), r, color=color)) + ax.set(xlabel=xlab, ylabel=ylab, aspect="equal") + ax.grid() + + return ax + + +def plot_optimal_layout( + boundaries=[], + x_final=[], + y_final=[], + x_init=[], + y_init=[], + D=126.0, + color_initial="tab:blue", + color_final="tab:orange", + norm=True, + ax=None, +): + """ + Plot the initial and final solution of a layout optimization study. + The turbine markers are properly scaled relative to the domain. + + Args: + boundaries (list(tuple)): boundary vertices in the form + [(x0,y0), (x1,y1), ... , (xN,yN)] + x_final (numpy.array(float)): x-positions of each turbine + in the optimal solution [m] + y_final (numpy.array(float)): y-positions of each turbine + in the optimal solution [m] + x_init (numpy.array(float)): x-positions of each turbine + in the initial solution [m] + y_init (numpy.array(float)): y-positions of each turbine + in the initial solution [m] + D (float): rotor diameter [m] + norm (bool): dictates whether the plot should be scaled by + rotor diameter. Defaults to True. + ax (:py:class:`matplotlib.pyplot.axes`, optional): axis to + plot layout + + Returns: + ax (:py:class:`matplotlib.pyplot.axes`, optional): axis + after plotting and formatting + + """ + + if ax is None: + _, ax = plt.subplots() + + if norm: + x1 = x_final / D + y1 = y_final / D + verts = np.array(boundaries) / D + r = 0.5 + xlab = "x/D" + ylab = "y/D" + + else: + x1 = x_final + y1 = y_final + verts = boundaries + r = D + xlab = "x [m]" + ylab = "y [m]" + + # Plot turbine locations + ax.scatter(x1, y1, s=0.01, color=color_final) + for x, y in zip(x1, y1): + ax.add_patch(plt.Circle((x, y), r, color=color_final)) + ax.set(xlabel=xlab, ylabel=ylab, aspect="equal") + ax.legend(["Final"], markerscale=50) + ax.grid() + + if verts.ndim == 2: + # Plot plant boundary for single boundary + for i in range(len(verts)): + if i == len(verts) - 1: + ax.plot([verts[i][0], verts[0][0]], [verts[i][1], verts[0][1]], "black") + else: + ax.plot( + [verts[i][0], verts[i + 1][0]], + [verts[i][1], verts[i + 1][1]], + "black", + ) + + else: + # Plot plant boundary for multiple boundaries + for ii in range(len(verts)): + sub_verts = verts[ii] + for i in range(len(sub_verts)): + if i == len(sub_verts) - 1: + ax.plot( + [sub_verts[i][0], sub_verts[0][0]], + [sub_verts[i][1], sub_verts[0][1]], + "black", + ) + else: + ax.plot( + [sub_verts[i][0], sub_verts[i + 1][0]], + [sub_verts[i][1], sub_verts[i + 1][1]], + "black", + ) + + return ax + + +########################################################################### +# Optimization performance methods +########################################################################### + + +def animate_layout_history( + filename=None, + layout_x=[], + layout_y=[], + boundaries=[], + D=126.0, + norm=True, + show=True, +): + """ + Animate the history of the wind farm layout. Plots the wind farm + layout at each iteration and saves to the given file as an MP4. + + Args: + filename (str, '.mp4'): name of animation file + layout_x (tuple(numpy.array)): container of x-positions at each + iteration. + layout_y (tuple(numpy.array)): container of y-positions at each + iteration. + boundaries (list(tuple)): boundary vertices in the form + [(x0,y0), (x1,y1), ... , (xN,yN)] + D (float): rotor diameter [m] + norm (bool): dictates whether the plot should be scaled by + rotor diameter. Defaults to True. + show (bool): dictates whether the animation is displayed before + closing. Defaults to True. + + """ + + if ax is None: + _, ax = plt.subplots() + + if filename is None: + raise ValueError("Must supply file name for layout history animation.") + + if norm: + x = layout_x / D + y = layout_y / D + verts = boundaries / D + xlab = "x/D" + ylab = "y/D" + + else: + x = layout_x + y = layout_y + verts = boundaries + xlab = "x [m]" + ylab = "y [m]" + + # Layout animation + fig, ax = plt.subplots() + ax.set(xlabel=xlab, ylabel=ylab, aspect="equal") + ax.grid() + + for i in range(len(verts)): + if i == len(verts) - 1: + ax.plot([verts[i][0], verts[0][0]], [verts[i][1], verts[0][1]], "black") + else: + ax.plot( + [verts[i][0], verts[i + 1][0]], [verts[i][1], verts[i + 1][1]], "black" + ) + + (line,) = ax.plot([], [], "o") + + # Function to update turbine positions + def animate(i): + line.set_data(x[i], y[i]) + ax.set_title(str(i)) + return (line,) + + # Animation + ani = animation.FuncAnimation(fig, animate, frames=len(x), repeat=False) + ani.save(filename) + if show: + plt.show() + else: + plt.close(fig) + + +def plot_convergence_history( + aep=[], + optimality=[], + feasibility=[], + ax_aep=None, + ax_opt=None, + ax_feas=None, +): + """ + Plots the convergence history of AEP (objective function), + optimality, and feasibility. Axes should be supplied for any + metrics that are to be plotted. + + Args: + aep (list(float)): AEP at each major iteration + optimality (list(float)): SNOPT optimality at each major iteration + feasibility (list(float)): SNOPT feasibility at each major iteration + ax_aep (:py:class:`matplotlib.pyplot.axes`, optional): axis to + plot AEP history + ax_opt (:py:class:`matplotlib.pyplot.axes`, optional): axis to + plot optimality history + ax_feas (:py:class:`matplotlib.pyplot.axes`, optional): axis to + plot feasibility history + + """ + + # Objective plot + if len(aep) > 0: + if ax_aep is None: + _, ax_aep = plt.subplots() + + ax_aep.plot([elem / 1e9 for elem in aep]) + ax_aep.set(xlabel="Iteration", ylabel="AEP [GWh]") + ax_aep.grid(True) + + # Optimality plot + if len(optimality) > 0: + if ax_opt is None: + _, ax_opt = plt.subplots() + + ax_opt.semilogy(optimality) + ax_opt.set(xlabel="Iteration", ylabel="Optimality [-]") + ax_opt.grid(True) + + # Feasibility plot + if len(feasibility) > 0: + if ax_feas is None: + _, ax_feas = plt.subplots() + + ax_feas.semilogy(feasibility) + ax_feas.set(xlabel="Iteration", ylabel="Feasibility [-]") + ax_feas.grid(True) + + +## Legacy functions + + +def plot_constraints(ax_boundary, ax_spacing, boundary_constraint, spacing_constraint): + """ + Plots the convergence history of the objective function and the wind farm + layout (optional) + + Args: + ax: matplotlib axis handle to plot AEP history. + obj (list(float)): A list of AEP at each major iteration. + layout (tuple(float)): A list of wind farm (x,y) layout at each major iteration. + boundaries (list(float)): A list of the boundary vertices in the form + [(x0,y0), (x1,y1), ... , (xN,yN)]. + D (float): rotor diameter + filename (str): name of .mp4 animation of layout progression. + """ + + # Boundary constraint plot + for n in range(len(boundary_constraint)): + ax_boundary.plot(boundary_constraint[n], alpha=0.3) + ax_boundary.set(xlabel="Iteration", ylabel="Boundary Constraint") + ax_boundary.grid() + + # Spacing constraint plot + ax_spacing.plot(spacing_constraint) + ax_spacing.set(xlabel="Iteration", ylabel="Spacing Constraint") + ax_spacing.grid() + + +def plot_flow_field(fi, ax, bounds, pts=200, cmin=2, cmax=10): + """ + Plots a filled contour map of the annually-averaged flow field. + + Args: + fi: FLOWERS interface with valid Fourier coefficients + ax: matplotlib axis handle to plot colormesh + bounds (np.array): domain limits in the form + ([x_min, x_max], [y_min, y_max]) + pts: grid resolution (uniform in x and y) + cmin: minimum wind speed for colorbar [m/s] + cmax: maximum wind speed for colorbar [m/s] + + """ + + # Enforce that fourier_coefficients() have been computed + if not hasattr(fi, "fs"): + print("Error, must compute Fourier coefficients before calculating wake") + return None + + # Define 2D grid based on defined domain limits + xx = np.linspace(bounds[0, 0], bounds[0, 1], pts) + yy = np.linspace(bounds[1, 0], bounds[1, 1], pts) + XX, YY = np.meshgrid(xx, yy) + + # Freestream velocity component + u0 = fi.fs.a_free[0] * np.pi + + # Superimpose wake velocity component from each turbine + for i in range(len(fi.layout_x)): + if i == 0: + du = fi.calculate_wake(XX - fi.layout_x[i], YY - fi.layout_y[i]) + else: + du += fi.calculate_wake(XX - fi.layout_x[i], YY - fi.layout_y[i]) + + # Compute average wake velocity + u = u0 - du + + # Mask points within rotor swept area + zz = np.logical_or( + np.isnan(u), + np.sqrt((XX - fi.layout_x[0]) ** 2 + (YY - fi.layout_y[0]) ** 2) < fi.D / 2, + ) + for j in range(len(fi.layout_x) - 1): + zz = np.logical_or( + zz, + np.sqrt((XX - fi.layout_x[j + 1]) ** 2 + (YY - fi.layout_y[j + 1]) ** 2) + < fi.D / 2, + ) + u_masked = np.ma.masked_where(zz, u) + + # Plot wake velocity colormesh and rotor swept areas + im = ax.pcolormesh(XX, YY, u_masked, vmin=cmin, vmax=cmax, cmap="coolwarm") + ax.plot(fi.layout_x, fi.layout_y, "ow") + + return im + + +def plot_floris_field(fli, ax, wind_rose, bounds, pts=200, cmin=2, cmax=10): + """ + Plots a filled contour map of the annually-averaged flow field using FLORIS. + + Args: + fli: FLORIS interface + fi: FLOWERS interface with valid Fourier coefficients + ax: matplotlib axis handle to plot colormesh + bounds (np.array): domain limits in the form + ([x_min, x_max], [y_min, y_max]) + pts: grid resolution (uniform in x and y) + cmin: minimum wind speed for colorbar [m/s] + cmax: maximum wind speed for colorbar [m/s] + """ + + # Resample wind rose by average wind speed to speed up plotting + wr = wind_rose.copy(deep=True) + wr = tl.resample_wind_direction(wr, wd=np.arange(0, 360, 5.0)) + wr = tl.resample_average_ws_by_wd(wr) + + # Numerically integrate over wind rose + for i in range(len(wr.wd)): # range(len(wr.wd)) + + # Redefine flow field for given wind speed and direction + fli.reinitialize(wind_directions=[wr.wd[i]], wind_speeds=[wr.ws[i]]) + if i == 0: + hor_plane = fli.calculate_horizontal_plane( + height=90.0, + x_resolution=pts, + y_resolution=pts, + x_bounds=[bounds[0, 0], bounds[0, 1]], + y_bounds=[bounds[1, 0], bounds[1, 1]], + ) + hor_plane.df.u = wr.freq_val[i] * hor_plane.df.u + else: + hor_plane1 = fli.calculate_horizontal_plane( + height=90.0, + x_resolution=pts, + y_resolution=pts, + x_bounds=[bounds[0, 0], bounds[0, 1]], + y_bounds=[bounds[1, 0], bounds[1, 1]], + ) + hor_plane.df.u += wr.freq_val[i] * hor_plane1.df.u + + # Reshape mesh grid for plotting + x1_mesh = hor_plane.df.x1.values.reshape( + hor_plane.resolution[1], hor_plane.resolution[0] + ) + x2_mesh = hor_plane.df.x2.values.reshape( + hor_plane.resolution[1], hor_plane.resolution[0] + ) + u_mesh = hor_plane.df.u.values.reshape( + hor_plane.resolution[1], hor_plane.resolution[0] + ) + + # Mask rotor swept areas (add rotor diameter variable) + zz = np.logical_or( + np.isnan(u_mesh), + np.sqrt( + (x1_mesh - fli.floris.farm.layout_x[0]) ** 2 + + (x2_mesh - fli.floris.farm.layout_y[0]) ** 2 + ) + < 126 / 2, + ) + for j in range(len(fli.floris.farm.layout_x) - 1): + zz = np.logical_or( + zz, + np.sqrt( + (x1_mesh - fli.floris.farm.layout_x[j + 1]) ** 2 + + (x2_mesh - fli.floris.farm.layout_y[j + 1]) ** 2 + ) + < 126 / 2, + ) + u = np.ma.masked_where(zz, u_mesh) + + # Plot wake velocity colormesh and rotor swept areas + im = ax.pcolormesh(x1_mesh, x2_mesh, u, cmap="coolwarm") # vmin=cmin, vmax=cmax, + for j in range(len(fli.floris.farm.layout_x)): + ax.scatter(fli.floris.farm.layout_x[j], fli.floris.farm.layout_y[j], c="white") + + return im diff --git a/pyproject.toml b/pyproject.toml index 024e5919..d174f24c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -11,7 +11,7 @@ build-backend = "setuptools.build_meta" # https://setuptools.pypa.io/en/latest/userguide/package_discovery.html [tool.setuptools.packages.find] where = ["."] -include = ["ard", "ard.*"] +include = ["ard", "ard.*", "flowers"] [tool.setuptools.package-data] "ard.api.default_systems" = ["*.yaml"] @@ -23,6 +23,7 @@ authors = [ {name = "Cory Frontin", email = "cory.frontin@nrel.gov"}, {name = "Rafael Mudafort", email = "rafael.mudafort@nrel.gov"}, {name = "Jared Thomas", email = "jared.thomas@nrel.gov"}, + {name = "Christopher Bay", email = "christopher.bay@nrel.gov"}, ] description = "A package for multidisciplinary and/or multifidelity wind farm design" readme= "README.md" @@ -39,6 +40,9 @@ classifiers = [ ] dependencies = [ "numpy", + "matplotlib", + "scipy", + "pandas", "floris>=4.3", "wisdem>=4.0.5", "NLopt", diff --git a/pyproject_flowers.toml b/pyproject_flowers.toml new file mode 100644 index 00000000..30108eed --- /dev/null +++ b/pyproject_flowers.toml @@ -0,0 +1,57 @@ + +# Project configuration for FLOWERS +# TOML spec: https://toml.io/en/v1.0.0 + +[build-system] +requires = [ + "setuptools", +] +build-backend = "setuptools.build_meta" + +# https://setuptools.pypa.io/en/latest/userguide/package_discovery.html +[tool.setuptools] +packages = ["flowers"] + +[project] +name = "FLOWERS" +version = "1.0.0-alpha0" +authors = [ +] +description = "A package for rapid energy simulation of wind farms" +readme= "README.md" +requires-python = ">=3.10, <3.13" +classifiers = [ + "Development Status :: 3 - Alpha", + "License :: OSI Approved :: Apache Software License", + "Operating System :: OS Independent", + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", + "Programming Language :: Python :: 3.12", + # "Programming Language :: Python :: 3.13", +] +dependencies = [ + "numpy", + "matplotlib", + "scipy", + "pandas", + "shapely", +] +[project.optional-dependencies] +dev = [ + "black", + "pytest", + "pytest-cov", +] +docs = [ + "pyxdsm", + "jupyter-book", + "sphinx-book-theme", + "sphinx-autodoc-typehints", +] +[project.urls] +# Homepage = "https://example.com" +Documentation = "https://wisdem.github.io/flowers" +Repository = "https://github.com/wisdem/flowers.git" +Issues = "https://github.com/wisdem/flowers/issues" +# Changelog = "https://github.com/wisdem/flowers/blob/main/CHANGELOG.md" \ No newline at end of file diff --git a/test/ard/system/api/inputs_offshore_floating/ard_system.yaml b/test/ard/system/api/inputs_offshore_floating/ard_system.yaml index b9f47aaa..5e40bde5 100644 --- a/test/ard/system/api/inputs_offshore_floating/ard_system.yaml +++ b/test/ard/system/api/inputs_offshore_floating/ard_system.yaml @@ -12,7 +12,7 @@ modeling_options: y: [5.0, -5.0, -5.0, 5.0] energy_resource: name: "OFL system test resource" - wind_resource: !include ../../../../../examples/data/windIO-plant_wind-resource_wrg-example.yaml + wind_resource: !include ../../../../../examples/ard/data/windIO-plant_wind-resource_wrg-example.yaml wind_farm: name: "OFL system test farm" layouts: @@ -31,7 +31,7 @@ modeling_options: 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ] - turbine: !include ../../../../../examples/data/windIO-plant_turbine_IEA-22MW-284m-RWT.yaml + turbine: !include ../../../../../examples/ard/data/windIO-plant_turbine_IEA-22MW-284m-RWT.yaml electrical_substations: - electrical_substation: coordinates: diff --git a/test/ard/system/api/inputs_offshore_monopile/ard_system.yaml b/test/ard/system/api/inputs_offshore_monopile/ard_system.yaml index 7ae376b7..e5415755 100644 --- a/test/ard/system/api/inputs_offshore_monopile/ard_system.yaml +++ b/test/ard/system/api/inputs_offshore_monopile/ard_system.yaml @@ -12,7 +12,7 @@ modeling_options: y: [5.0, -5.0, -5.0, 5.0] energy_resource: name: "OFB system test resource" - wind_resource: !include ../../../../../examples/data/windIO-plant_wind-resource_wrg-example.yaml + wind_resource: !include ../../../../../examples/ard/data/windIO-plant_wind-resource_wrg-example.yaml wind_farm: name: "OFB system test farm" layouts: @@ -31,7 +31,7 @@ modeling_options: 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ] - turbine: !include ../../../../../examples/data/windIO-plant_turbine_IEA-22MW-284m-RWT.yaml + turbine: !include ../../../../../examples/ard/data/windIO-plant_turbine_IEA-22MW-284m-RWT.yaml electrical_substations: - electrical_substation: coordinates: diff --git a/test/ard/system/api/inputs_onshore/ard_system.yaml b/test/ard/system/api/inputs_onshore/ard_system.yaml index 517560e5..3c3202d4 100644 --- a/test/ard/system/api/inputs_onshore/ard_system.yaml +++ b/test/ard/system/api/inputs_onshore/ard_system.yaml @@ -12,7 +12,7 @@ modeling_options: y: [5.0, -5.0, -5.0, 5.0] energy_resource: name: "LB system test resource" - wind_resource: !include ../../../../../examples/data/windIO-plant_wind-resource_wrg-example.yaml + wind_resource: !include ../../../../../examples/ard/data/windIO-plant_wind-resource_wrg-example.yaml wind_farm: name: "LB system test farm" layouts: @@ -31,7 +31,7 @@ modeling_options: 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ] - turbine: !include ../../../../../examples/data/windIO-plant_turbine_IEA-3.4MW-130m-RWT.yaml + turbine: !include ../../../../../examples/ard/data/windIO-plant_turbine_IEA-3.4MW-130m-RWT.yaml electrical_substations: - electrical_substation: coordinates: diff --git a/test/ard/system/collection/test_optiwindnet.py b/test/ard/system/collection/test_optiwindnet.py index fe5ede69..03655ad6 100644 --- a/test/ard/system/collection/test_optiwindnet.py +++ b/test/ard/system/collection/test_optiwindnet.py @@ -29,12 +29,14 @@ def setup_method(self): filename_turbine = ( Path(ard.__file__).parents[1] / "examples" + / "ard" / "data" / "windIO-plant_turbine_IEA-3.4MW-130m-RWT.yaml" ) filename_windresource = ( Path(ard.__file__).parents[1] / "examples" + / "ard" / "data" / "windIO-plant_wind-resource_wrg-example.yaml" ) diff --git a/test/ard/system/cost/test_spacing_approximations_connections.py b/test/ard/system/cost/test_spacing_approximations_connections.py index 39df136c..ccd64bde 100644 --- a/test/ard/system/cost/test_spacing_approximations_connections.py +++ b/test/ard/system/cost/test_spacing_approximations_connections.py @@ -13,6 +13,7 @@ def setup_method(self): filename_turbine = ( Path(ard.__file__).parents[1] / "examples" + / "ard" / "data" / "windIO-plant_turbine_IEA-3.4MW-130m-RWT.yaml" ) diff --git a/test/ard/system/geometry/test_constraints.py b/test/ard/system/geometry/test_constraints.py index 0564ee3d..a43fc749 100644 --- a/test/ard/system/geometry/test_constraints.py +++ b/test/ard/system/geometry/test_constraints.py @@ -26,6 +26,7 @@ def setup_method(self): path_turbine = ( Path(ard.__file__).parents[1] / "examples" + / "ard" / "data" / "windIO-plant_turbine_IEA-3.4MW-130m-RWT.yaml" ) diff --git a/test/ard/system/offshore/test_mooring_packing.py b/test/ard/system/offshore/test_mooring_packing.py index 7e0b6bc0..fadba154 100644 --- a/test/ard/system/offshore/test_mooring_packing.py +++ b/test/ard/system/offshore/test_mooring_packing.py @@ -21,6 +21,7 @@ def setup_method(self): filename_turbine = ( Path(ard.__file__).parents[1] / "examples" + / "ard" / "data" / "windIO-plant_turbine_IEA-22MW-284m-RWT.yaml" ) # windIO turbine specification diff --git a/test/ard/unit/api/inputs_onshore/ard_system_bad_windio.yaml b/test/ard/unit/api/inputs_onshore/ard_system_bad_windio.yaml index e1b5dc5e..4a821de8 100644 --- a/test/ard/unit/api/inputs_onshore/ard_system_bad_windio.yaml +++ b/test/ard/unit/api/inputs_onshore/ard_system_bad_windio.yaml @@ -10,7 +10,7 @@ modeling_options: - x: [5.0, 5.0, -5.0, -5.0] energy_resource: name: "LB system test resource" - wind_resource: !include ../../../../../examples/data/windIO-plant_wind-resource_wrg-example.yaml + wind_resource: !include ../../../../../examples/ard/data/windIO-plant_wind-resource_wrg-example.yaml wind_farm: name: "LB system test farm" layouts: @@ -29,7 +29,7 @@ modeling_options: 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ] - turbine: !include ../../../../../examples/data/windIO-plant_turbine_IEA-3.4MW-130m-RWT.yaml + turbine: !include ../../../../../examples/ard/data/windIO-plant_turbine_IEA-3.4MW-130m-RWT.yaml electrical_substations: - electrical_substation: coordinates: diff --git a/test/ard/unit/api/inputs_onshore/windio.yaml b/test/ard/unit/api/inputs_onshore/windio.yaml index 91041a19..43dea33a 100644 --- a/test/ard/unit/api/inputs_onshore/windio.yaml +++ b/test/ard/unit/api/inputs_onshore/windio.yaml @@ -7,7 +7,7 @@ site: y: [ 3000.0, 1500.0, -1500.0, -3000.0, -3000.0, -1500.0, 1500.0, 3000.0] energy_resource: name: Ard Example 01 offshore energy resource - wind_resource: !include ../../../../../examples/data/windIO-plant_wind-resource_wrg-example.yaml + wind_resource: !include ../../../../../examples/ard/data/windIO-plant_wind-resource_wrg-example.yaml wind_farm: name: Ard Example 01 offshore wind farm layouts: @@ -26,7 +26,7 @@ wind_farm: 1250.0, 1250.0, 1250.0, 1250.0, 1250.0, 2500.0, 2500.0, 2500.0, 2500.0, 2500.0 ] - turbine: !include ../../../../../examples/data/windIO-plant_turbine_IEA-3.4MW-130m-RWT.yaml + turbine: !include ../../../../../examples/ard/data/windIO-plant_turbine_IEA-3.4MW-130m-RWT.yaml electrical_substations: - electrical_substation: coordinates: diff --git a/test/ard/unit/cost/test_orbit_wrap.py b/test/ard/unit/cost/test_orbit_wrap.py index 2e6e8766..d1d501c8 100644 --- a/test/ard/unit/cost/test_orbit_wrap.py +++ b/test/ard/unit/cost/test_orbit_wrap.py @@ -20,6 +20,7 @@ def test_raise_error(self): filename_turbine = ( Path(ard.__file__).parents[1] / "examples" + / "ard" / "data" / "windIO-plant_turbine_IEA-22MW-284m-RWT.yaml" ).absolute() # toolset generalized turbine specification @@ -202,6 +203,7 @@ def test_baseline_farm(self, subtests): filename_turbine = ( Path(ard.__file__).parents[1] / "examples" + / "ard" / "data" / "windIO-plant_turbine_IEA-22MW-284m-RWT.yaml" ).absolute() # toolset generalized turbine specification @@ -399,6 +401,7 @@ def setup_method(self): filename_turbine = ( Path(ard.__file__).parents[1] / "examples" + / "ard" / "data" / "windIO-plant_turbine_IEA-22MW-284m-RWT.yaml" ).absolute() # toolset generalized turbine specification diff --git a/test/ard/unit/cost/test_wisdem_wrap.py b/test/ard/unit/cost/test_wisdem_wrap.py index 8d4ae901..e200170b 100644 --- a/test/ard/unit/cost/test_wisdem_wrap.py +++ b/test/ard/unit/cost/test_wisdem_wrap.py @@ -18,6 +18,7 @@ def setup_method(self): filename_turbine = ( Path(ard.__file__).parents[1] / "examples" + / "ard" / "data" / "windIO-plant_turbine_IEA-3.4MW-130m-RWT.yaml" ) @@ -112,6 +113,7 @@ def setup_method(self): filename_turbine = ( Path(ard.__file__).parents[1] / "examples" + / "ard" / "data" / "windIO-plant_turbine_IEA-22MW-284m-RWT.yaml" ) diff --git a/test/ard/unit/farm_aero/test_floris.py b/test/ard/unit/farm_aero/test_floris.py index 387a7113..a956f5e0 100644 --- a/test/ard/unit/farm_aero/test_floris.py +++ b/test/ard/unit/farm_aero/test_floris.py @@ -41,6 +41,7 @@ def setup_method(self): path_turbine = ( Path(ard.__file__).parents[1] / "examples" + / "ard" / "data" / "windIO-plant_turbine_IEA-3.4MW-130m-RWT.yaml" ) @@ -174,6 +175,7 @@ def setup_method(self): path_turbine = ( Path(ard.__file__).parents[1] / "examples" + / "ard" / "data" / "windIO-plant_turbine_IEA-3.4MW-130m-RWT.yaml" ) diff --git a/test/ard/unit/geographic/test_geomorphology.py b/test/ard/unit/geographic/test_geomorphology.py index 05f9b35a..1b2d260e 100644 --- a/test/ard/unit/geographic/test_geomorphology.py +++ b/test/ard/unit/geographic/test_geomorphology.py @@ -335,6 +335,7 @@ def test_load_moorpy_bathymetry(self, subtests): file_bathy = ( Path(ard.__file__).parents[1] / "examples" + / "ard" / "data" / "offshore" / "GulfOfMaine_bathymetry_100x99.txt" diff --git a/test/flowers/__init__.py b/test/flowers/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/test/flowers/_archive/conftest.py b/test/flowers/_archive/conftest.py new file mode 100644 index 00000000..e69de29b diff --git a/test/flowers/_archive/regression/test_flowers_aep.py b/test/flowers/_archive/regression/test_flowers_aep.py new file mode 100644 index 00000000..09d1ac22 --- /dev/null +++ b/test/flowers/_archive/regression/test_flowers_aep.py @@ -0,0 +1,54 @@ +from pathlib import Path + +import numpy as np +import pandas as pd +from pytest import approx + +import flowers +from flowers import FlowersModel + + +def test_regression_aep(): + # Create generic layout + D = 126.0 + layout_x = D * np.array( + [ + 0.0, + 0.0, + 0.0, + 7.0, + 7.0, + 7.0, + 14.0, + 14.0, + 14.0, + 21.0, + 21.0, + 21.0, + 28.0, + 28.0, + 28.0, + ] + ) + layout_y = D * np.array( + [0.0, 7.0, 14.0, 0.0, 7.0, 14.0, 0.0, 7.0, 14.0, 0.0, 7.0, 14.0, 0.0, 7.0, 14.0] + ) + + # Load in wind data + wind_rose_file = Path( + Path(flowers.__file__).parent, + "..", + "examples", + "flowers", + "data", + "HKW_wind_rose.csv", + ) + df = pd.read_csv(wind_rose_file) + + # Setup FLOWERS model + flowers_model = FlowersModel(df, layout_x, layout_y) + + # Calculate the AEP + aep = flowers_model.calculate_aep() + + assert aep == approx(350958996034.8524, 1e-3) diff --git a/test/flowers/system/rotation/__init__.py b/test/flowers/system/rotation/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/test/flowers/system/rotation/precooked.yaml b/test/flowers/system/rotation/precooked.yaml new file mode 100644 index 00000000..9fe3bb83 --- /dev/null +++ b/test/flowers/system/rotation/precooked.yaml @@ -0,0 +1,1177 @@ +logging: + console: + enable: true + level: WARNING + file: + enable: false + level: WARNING +solver: + type: turbine_grid + turbine_grid_points: 3 +wake: + model_strings: + combination_model: sosfs + deflection_model: gauss + turbulence_model: crespo_hernandez + velocity_model: gauss + enable_secondary_steering: true + enable_yaw_added_recovery: true + enable_active_wake_mixing: false + enable_transverse_velocities: true + wake_deflection_parameters: + gauss: + ad: 0.0 + alpha: 0.58 + bd: 0.0 + beta: 0.077 + dm: 1.0 + ka: 0.38 + kb: 0.004 + jimenez: + ad: 0.0 + bd: 0.0 + kd: 0.05 + empirical_gauss: + horizontal_deflection_gain_D: 3.0 + vertical_deflection_gain_D: -1 + deflection_rate: 22 + mixing_gain_deflection: 0.0 + yaw_added_mixing_gain: 0.0 + wake_turbulence_parameters: + crespo_hernandez: + initial: 0.1 + constant: 0.5 + ai: 0.8 + downstream: -0.32 + wake_induced_mixing: + atmospheric_ti_gain: 0.0 + wake_velocity_parameters: + gauss: + alpha: 0.58 + beta: 0.077 + ka: 0.38 + kb: 0.004 + jensen: + we: 0.05 + cc: + a_s: 0.179367259 + b_s: 0.0118889215 + c_s1: 0.0563691592 + c_s2: 0.13290157 + a_f: 3.11 + b_f: -0.68 + c_f: 2.41 + alpha_mod: 1.0 + turbopark: + A: 0.04 + sigma_max_rel: 4.0 + turboparkgauss: + A: 0.04 + include_mirror_wake: true + empirical_gauss: + wake_expansion_rates: + - 0.023 + - 0.008 + breakpoints_D: + - 10 + sigma_0_D: 0.28 + smoothing_length_D: 2.0 + mixing_gain_velocity: 2.0 + awc_wake_exp: 1.2 + awc_wake_denominator: 400 +farm: + layout_x: + - 159.4005622195224 + - 799.839360501436 + - 510.5330763394482 + layout_y: + - -987.8891420375654 + - -344.73829003042766 + - -1067.7019703118387 + turbine_type: + - turbine_type: IEA-3.4MW-130m-RWT + hub_height: 110.0 + rotor_diameter: 130.0 + TSR: 8.01754386 + operation_model: cosine-loss + power_thrust_table: + ref_air_density: 1.225 + ref_tilt: 5.0 + cosine_loss_exponent_yaw: 1.88 + cosine_loss_exponent_tilt: 1.88 + wind_speed: + - 0.0 + - 2.9999 + - 3.0 + - 3.539159827293725 + - 4.048452625194753 + - 4.526747636329809 + - 4.9729829261640415 + - 5.386167740761253 + - 5.76538470650541 + - 6.109791866899474 + - 6.418624551919384 + - 6.691197075772739 + - 6.9269042592927175 + - 7.125222773587183 + - 7.285712301959686 + - 7.408016517522657 + - 7.491863874332182 + - 7.537068210287934 + - 7.543529160459564 + - 7.575825940997163 + - 7.646808745314545 + - 7.756319973835029 + - 7.9041164842451614 + - 8.089870131331459 + - 8.313168495544858 + - 8.57351579867522 + - 8.870334004602977 + - 9.202964102683834 + - 9.570667570917234 + - 9.812675420388173 + - 10.407952984173718 + - 10.875675946194342 + - 11.374758439769723 + - 11.904092376955541 + - 12.462502504037497 + - 13.048749010888468 + - 13.661530283657028 + - 14.299485794675704 + - 14.96119912317263 + - 15.645201100079884 + - 16.349973069956214 + - 17.0739502627819 + - 17.815525268139492 + - 18.57305160406695 + - 19.344847372659316 + - 20.12919899430266 + - 20.924365012249385 + - 21.728579959087647 + - 22.54005827652058 + - 23.356998279752162 + - 24.177586157678096 + - 25.0 + - 25.01 + - 50.0 + power: + - 0.0 + - 0.0 + - 55.485650388273946 + - 132.60594949759707 + - 229.20900366592642 + - 342.712543238107 + - 468.6008480810115 + - 602.8982926864225 + - 740.4915214111589 + - 880.4618007632126 + - 1020.8388498653875 + - 1156.4922707088438 + - 1283.0654945862623 + - 1396.4539297233146 + - 1492.9570847061225 + - 1569.4124611097482 + - 1623.3079057292164 + - 1652.8696585473417 + - 1657.1239382759463 + - 1678.4995701413168 + - 1726.1238569674936 + - 1801.3513158108537 + - 1906.3000805443926 + - 2043.8824757534887 + - 2217.8441255643593 + - 2432.8096396619285 + - 2694.3325007912877 + - 3008.9464213493393 + - 3384.21515097733 + - 3622.320140431238 + - 3622.3451711827993 + - 3622.3459040961607 + - 3622.3452257887666 + - 3622.3451528700944 + - 3622.3574440121315 + - 3622.345165812334 + - 3622.3469910649405 + - 3622.3544214545163 + - 3622.3813238557404 + - 3622.3451993631197 + - 3622.3453899482156 + - 3622.3467481326475 + - 3622.3475934972626 + - 3622.347589398676 + - 3622.3647319812258 + - 3622.3816921026446 + - 3622.350066618316 + - 3622.3450868531027 + - 3622.345851188797 + - 3622.349377884576 + - 3622.3452304136345 + - 3622.457915370389 + - 0.0 + - 0.0 + thrust_coefficient: + - 0.0 + - 0.0 + - 0.8140283854894109 + - 0.8041771483638767 + - 0.7982113770373608 + - 0.7945189424815273 + - 0.791309102845844 + - 0.786112395855186 + - 0.7760330855568125 + - 0.766405559051694 + - 0.766405559051694 + - 0.7664055590516942 + - 0.7664055590516942 + - 0.766405559051694 + - 0.766405559051694 + - 0.7664055590516939 + - 0.7664055590516942 + - 0.766405559051694 + - 0.7664055590516943 + - 0.766405559051694 + - 0.7664055590516939 + - 0.7664055590516939 + - 0.766405559051694 + - 0.7664055590516943 + - 0.7664055590516939 + - 0.7664055590516944 + - 0.766405559051694 + - 0.7664055590516939 + - 0.7664055590516942 + - 0.8067632947740724 + - 0.5306091101592825 + - 0.4466865448824549 + - 0.3801106706360693 + - 0.3254312907353938 + - 0.2798429968252726 + - 0.2414778385628342 + - 0.20904693927376936 + - 0.1815926995732083 + - 0.1582789238889643 + - 0.13843565625820317 + - 0.12141358275643703 + - 0.10691750603387402 + - 0.09454669711591918 + - 0.08392954309933542 + - 0.07479548070995087 + - 0.06691762496259271 + - 0.060104597639619645 + - 0.05418408689005904 + - 0.04906125901748401 + - 0.044599887927284615 + - 0.04063092405931261 + - 0.03722282641218611 + - 0.0 + - 0.0 +flow_field: + wind_speeds: + - 1.0 + - 2.0 + - 3.0 + - 4.0 + - 5.0 + - 6.0 + - 7.0 + - 8.0 + - 9.0 + - 10.0 + - 11.0 + - 12.0 + - 13.0 + - 14.0 + - 15.0 + - 16.0 + - 17.0 + - 18.0 + - 19.0 + - 20.0 + - 21.0 + - 22.0 + - 23.0 + - 24.0 + - 25.0 + - 1.0 + - 2.0 + - 3.0 + - 4.0 + - 5.0 + - 6.0 + - 7.0 + - 8.0 + - 9.0 + - 10.0 + - 11.0 + - 12.0 + - 13.0 + - 14.0 + - 15.0 + - 16.0 + - 17.0 + - 18.0 + - 19.0 + - 20.0 + - 21.0 + - 22.0 + - 23.0 + - 24.0 + - 25.0 + - 1.0 + - 2.0 + - 3.0 + - 4.0 + - 5.0 + - 6.0 + - 7.0 + - 8.0 + - 9.0 + - 10.0 + - 11.0 + - 12.0 + - 13.0 + - 14.0 + - 15.0 + - 16.0 + - 17.0 + - 18.0 + - 19.0 + - 20.0 + - 21.0 + - 22.0 + - 23.0 + - 24.0 + - 25.0 + - 1.0 + - 2.0 + - 3.0 + - 4.0 + - 5.0 + - 6.0 + - 7.0 + - 8.0 + - 9.0 + - 10.0 + - 11.0 + - 12.0 + - 13.0 + - 14.0 + - 15.0 + - 16.0 + - 17.0 + - 18.0 + - 19.0 + - 20.0 + - 21.0 + - 22.0 + - 23.0 + - 24.0 + - 25.0 + - 1.0 + - 2.0 + - 3.0 + - 4.0 + - 5.0 + - 6.0 + - 7.0 + - 8.0 + - 9.0 + - 10.0 + - 11.0 + - 12.0 + - 13.0 + - 14.0 + - 15.0 + - 16.0 + - 17.0 + - 18.0 + - 19.0 + - 20.0 + - 21.0 + - 22.0 + - 23.0 + - 24.0 + - 25.0 + - 1.0 + - 2.0 + - 3.0 + - 4.0 + - 5.0 + - 6.0 + - 7.0 + - 8.0 + - 9.0 + - 10.0 + - 11.0 + - 12.0 + - 13.0 + - 14.0 + - 15.0 + - 16.0 + - 17.0 + - 18.0 + - 19.0 + - 20.0 + - 21.0 + - 22.0 + - 23.0 + - 24.0 + - 25.0 + - 1.0 + - 2.0 + - 3.0 + - 4.0 + - 5.0 + - 6.0 + - 7.0 + - 8.0 + - 9.0 + - 10.0 + - 11.0 + - 12.0 + - 13.0 + - 14.0 + - 15.0 + - 16.0 + - 17.0 + - 18.0 + - 19.0 + - 20.0 + - 21.0 + - 22.0 + - 23.0 + - 24.0 + - 25.0 + - 1.0 + - 2.0 + - 3.0 + - 4.0 + - 5.0 + - 6.0 + - 7.0 + - 8.0 + - 9.0 + - 10.0 + - 11.0 + - 12.0 + - 13.0 + - 14.0 + - 15.0 + - 16.0 + - 17.0 + - 18.0 + - 19.0 + - 20.0 + - 21.0 + - 22.0 + - 23.0 + - 24.0 + - 25.0 + - 1.0 + - 2.0 + - 3.0 + - 4.0 + - 5.0 + - 6.0 + - 7.0 + - 8.0 + - 9.0 + - 10.0 + - 11.0 + - 12.0 + - 13.0 + - 14.0 + - 15.0 + - 16.0 + - 17.0 + - 18.0 + - 19.0 + - 20.0 + - 21.0 + - 22.0 + - 23.0 + - 24.0 + - 25.0 + - 1.0 + - 2.0 + - 3.0 + - 4.0 + - 5.0 + - 6.0 + - 7.0 + - 8.0 + - 9.0 + - 10.0 + - 11.0 + - 12.0 + - 13.0 + - 14.0 + - 15.0 + - 16.0 + - 17.0 + - 18.0 + - 19.0 + - 20.0 + - 21.0 + - 22.0 + - 23.0 + - 24.0 + - 25.0 + - 1.0 + - 2.0 + - 3.0 + - 4.0 + - 5.0 + - 6.0 + - 7.0 + - 8.0 + - 9.0 + - 10.0 + - 11.0 + - 12.0 + - 13.0 + - 14.0 + - 15.0 + - 16.0 + - 17.0 + - 18.0 + - 19.0 + - 20.0 + - 21.0 + - 22.0 + - 23.0 + - 24.0 + - 25.0 + - 1.0 + - 2.0 + - 3.0 + - 4.0 + - 5.0 + - 6.0 + - 7.0 + - 8.0 + - 9.0 + - 10.0 + - 11.0 + - 12.0 + - 13.0 + - 14.0 + - 15.0 + - 16.0 + - 17.0 + - 18.0 + - 19.0 + - 20.0 + - 21.0 + - 22.0 + - 23.0 + - 24.0 + - 25.0 + wind_directions: + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - 30.0 + - 30.0 + - 30.0 + - 30.0 + - 30.0 + - 30.0 + - 30.0 + - 30.0 + - 30.0 + - 30.0 + - 30.0 + - 30.0 + - 30.0 + - 30.0 + - 30.0 + - 30.0 + - 30.0 + - 30.0 + - 30.0 + - 30.0 + - 30.0 + - 30.0 + - 30.0 + - 30.0 + - 30.0 + - 60.0 + - 60.0 + - 60.0 + - 60.0 + - 60.0 + - 60.0 + - 60.0 + - 60.0 + - 60.0 + - 60.0 + - 60.0 + - 60.0 + - 60.0 + - 60.0 + - 60.0 + - 60.0 + - 60.0 + - 60.0 + - 60.0 + - 60.0 + - 60.0 + - 60.0 + - 60.0 + - 60.0 + - 60.0 + - 90.0 + - 90.0 + - 90.0 + - 90.0 + - 90.0 + - 90.0 + - 90.0 + - 90.0 + - 90.0 + - 90.0 + - 90.0 + - 90.0 + - 90.0 + - 90.0 + - 90.0 + - 90.0 + - 90.0 + - 90.0 + - 90.0 + - 90.0 + - 90.0 + - 90.0 + - 90.0 + - 90.0 + - 90.0 + - 120.0 + - 120.0 + - 120.0 + - 120.0 + - 120.0 + - 120.0 + - 120.0 + - 120.0 + - 120.0 + - 120.0 + - 120.0 + - 120.0 + - 120.0 + - 120.0 + - 120.0 + - 120.0 + - 120.0 + - 120.0 + - 120.0 + - 120.0 + - 120.0 + - 120.0 + - 120.0 + - 120.0 + - 120.0 + - 150.0 + - 150.0 + - 150.0 + - 150.0 + - 150.0 + - 150.0 + - 150.0 + - 150.0 + - 150.0 + - 150.0 + - 150.0 + - 150.0 + - 150.0 + - 150.0 + - 150.0 + - 150.0 + - 150.0 + - 150.0 + - 150.0 + - 150.0 + - 150.0 + - 150.0 + - 150.0 + - 150.0 + - 150.0 + - 180.0 + - 180.0 + - 180.0 + - 180.0 + - 180.0 + - 180.0 + - 180.0 + - 180.0 + - 180.0 + - 180.0 + - 180.0 + - 180.0 + - 180.0 + - 180.0 + - 180.0 + - 180.0 + - 180.0 + - 180.0 + - 180.0 + - 180.0 + - 180.0 + - 180.0 + - 180.0 + - 180.0 + - 180.0 + - 210.0 + - 210.0 + - 210.0 + - 210.0 + - 210.0 + - 210.0 + - 210.0 + - 210.0 + - 210.0 + - 210.0 + - 210.0 + - 210.0 + - 210.0 + - 210.0 + - 210.0 + - 210.0 + - 210.0 + - 210.0 + - 210.0 + - 210.0 + - 210.0 + - 210.0 + - 210.0 + - 210.0 + - 210.0 + - 240.0 + - 240.0 + - 240.0 + - 240.0 + - 240.0 + - 240.0 + - 240.0 + - 240.0 + - 240.0 + - 240.0 + - 240.0 + - 240.0 + - 240.0 + - 240.0 + - 240.0 + - 240.0 + - 240.0 + - 240.0 + - 240.0 + - 240.0 + - 240.0 + - 240.0 + - 240.0 + - 240.0 + - 240.0 + - 270.0 + - 270.0 + - 270.0 + - 270.0 + - 270.0 + - 270.0 + - 270.0 + - 270.0 + - 270.0 + - 270.0 + - 270.0 + - 270.0 + - 270.0 + - 270.0 + - 270.0 + - 270.0 + - 270.0 + - 270.0 + - 270.0 + - 270.0 + - 270.0 + - 270.0 + - 270.0 + - 270.0 + - 270.0 + - 300.0 + - 300.0 + - 300.0 + - 300.0 + - 300.0 + - 300.0 + - 300.0 + - 300.0 + - 300.0 + - 300.0 + - 300.0 + - 300.0 + - 300.0 + - 300.0 + - 300.0 + - 300.0 + - 300.0 + - 300.0 + - 300.0 + - 300.0 + - 300.0 + - 300.0 + - 300.0 + - 300.0 + - 300.0 + - 330.0 + - 330.0 + - 330.0 + - 330.0 + - 330.0 + - 330.0 + - 330.0 + - 330.0 + - 330.0 + - 330.0 + - 330.0 + - 330.0 + - 330.0 + - 330.0 + - 330.0 + - 330.0 + - 330.0 + - 330.0 + - 330.0 + - 330.0 + - 330.0 + - 330.0 + - 330.0 + - 330.0 + - 330.0 + wind_veer: 0.0 + wind_shear: 0.12 + air_density: 1.225 + turbulence_intensities: + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + - 0.1 + reference_wind_height: 90.0 +name: GCH +description: 'Default initialization: Gauss-Curl hybrid model (GCH)' +floris_version: v4 diff --git a/test/flowers/system/rotation/rotational_consistency.npz b/test/flowers/system/rotation/rotational_consistency.npz new file mode 100644 index 00000000..02a0505b Binary files /dev/null and b/test/flowers/system/rotation/rotational_consistency.npz differ diff --git a/test/flowers/system/rotation/sandbox.py b/test/flowers/system/rotation/sandbox.py new file mode 100644 index 00000000..ecf64ce8 --- /dev/null +++ b/test/flowers/system/rotation/sandbox.py @@ -0,0 +1,194 @@ +#!/usr/bin/env python + +# import and set up packages +from pathlib import Path + +import tqdm + +import numpy as np +import matplotlib.pyplot as plt + +import ard +import floris + +import test_rotational_workbench + +plt.style.use(ard.get_house_style(use_tex=True)) + +## configuration settings + +path_floris_config = Path.cwd() / "precooked.yaml" +show_plots = False +save_plots = True +save_data = True +run_debug_mode = False + +## generate demo plots for the testing workbench setup + +if show_plots or save_plots: + # demo wraparound gaussian pulse resource for nearly single-directional wind + (WD, WS, FREQ), df_wr = test_rotational_workbench.f_windrose(-5.0, 8.0) + fig, ax = plt.subplots() + ct0 = ax.contourf(WD, WS, FREQ) + ax.set_xlabel("wind direction, (°)") + ax.set_ylabel("wind speed, (m/s)") + cb = fig.colorbar(ct0) + cb.set_label("probability density (-)") + if save_plots: + fig.savefig("wind_resource_demo.png", dpi=300, bbox_inches="tight") + plt.show() if show_plots else plt.close(fig) + + # demo the layout generator w/o rotational symmetries + fig, ax = plt.subplots() + ax.scatter( + *[ + v.flat + for v in test_rotational_workbench.layout_generator( + diameter=130.0, orientation=0.0 + ) + ] + ) + ax.set_xlabel("turbine x location/relative easting, (m)") + ax.set_ylabel("turbine y location/relative northing, (m)") + if save_plots: + fig.savefig("layout_demo.png", dpi=300, bbox_inches="tight") + plt.show() if show_plots else plt.close(fig) + +## set up the aerodynamic models and call functions + +# generate the FLORIS model w/ the pre-cooked IEA 3.4 RWT model +floris_model = floris.FlorisModel(path_floris_config) + + +# a function to run either wind direction or orientation sweeps +def run_floris(wd_val=0.0, orientation=0.0): + WD, WS, FREQ = test_rotational_workbench.f_windrose( + wd_center=wd_val, + ws_center=8.0, + )[0] + wind_rose = floris.wind_data.WindRose( + wind_directions=WD[:, 0], + wind_speeds=WS[0, :], + ti_table=0.06, + freq_table=FREQ, + ) + + X, Y = test_rotational_workbench.layout_generator(D_rotor, orientation) + floris_model.set( + wind_data=wind_rose, + layout_x=X.flatten(), + layout_y=Y.flatten(), + ) + floris_model.run() + return floris_model.get_farm_AEP() / 1e9 # report in GWh + + +# extract FLORIS data for FLOWERS s.t. the two models are 1:1 +turbine_def = floris_model.core.farm.turbine_definitions[0] +D_rotor = turbine_def["rotor_diameter"] +# thrust and power curves +u_cp = u_ct = np.array(turbine_def["power_thrust_table"]["wind_speed"][1:]) +cp = ( + np.array(turbine_def["power_thrust_table"]["power"][1:]) + * 1e3 + / ( + 0.5 + * turbine_def["power_thrust_table"]["ref_air_density"] + * (np.pi * D_rotor**2 / 4) + * u_cp**3 + ) +) +ct = np.array(turbine_def["power_thrust_table"]["thrust_coefficient"][1:]) +# detect cutin and cutout +U_cutin = float( + u_cp[ + np.squeeze( + np.argwhere( + np.isclose(cp, 0.0) & (np.array(list(range(len(cp)))) < len(cp) / 2) + )[0] + ) + + 1 + ] +) +U_cutout = float( + u_cp[ + np.squeeze( + np.argwhere( + np.isclose(cp, 0.0) & (np.array(list(range(len(cp)))) > len(cp) / 2) + )[0] + ) + - 1 + ] +) +print(f"detected cut in {U_cutin} and cut out {U_cutout}. verify.") + +# pack up flowers model +flowers_turbine = { + "D": D_rotor, + "cp": cp, + "u_cp": u_cp, + "ct": ct, + "u_ct": u_ct, + "U": U_cutout, +} + +## run an experiment: move the pulse wind resource around the compass + +# generate and run the cases +wd_vec = np.arange(0.0, 360.001, 1.0 if not run_debug_mode else 10.0) +AEP_FLOWERS_e1_vec = np.zeros_like(wd_vec) +AEP_FLORIS_e1_vec = np.zeros_like(wd_vec) +for idx, wd_val in enumerate(tqdm.tqdm(wd_vec)): + AEP_FLOWERS_e1_vec[idx] = test_rotational_workbench.run_FLOWERS( + flowers_turbine, wd_val=wd_val + ) + AEP_FLORIS_e1_vec[idx] = run_floris(wd_val=wd_val) + +# plot the results of the experiment +if show_plots or save_plots: + fig, axes = plt.subplots(2, 1) + axes[0].plot(wd_vec, AEP_FLOWERS_e1_vec, label="FLOWERS") + axes[0].plot(wd_vec, AEP_FLORIS_e1_vec, label="FLORIS") + axes[0].legend() + axes[1].hist((AEP_FLOWERS_e1_vec - AEP_FLORIS_e1_vec) / AEP_FLORIS_e1_vec * 100) + axes[1].set_xlabel("percent difference in FLOWERS w.r.t. FLORIS, (\\%)") + if save_plots: + fig.savefig("compare_AEP_vs_WD.png", dpi=300, bbox_inches="tight") + plt.show() if show_plots else plt.close(fig) + +## run an experiment: rotate the farm around the compass w/ fixed resource + +# generate and run the cases +orientation_vec = np.arange(0.0, 360.001, 1.0 if not run_debug_mode else 10.0) +AEP_FLOWERS_e2_vec = np.zeros_like(orientation_vec) +AEP_FLORIS_e2_vec = np.zeros_like(orientation_vec) +for idx, orientation_val in enumerate(tqdm.tqdm(orientation_vec)): + AEP_FLOWERS_e2_vec[idx] = test_rotational_workbench.run_FLOWERS( + flowers_turbine, orientation=orientation_val + ) + AEP_FLORIS_e2_vec[idx] = run_floris(orientation=orientation_val) + +# plot the results of the experiment +if show_plots or save_plots: + fig, axes = plt.subplots(2, 1) + axes[0].plot(orientation_vec, AEP_FLOWERS_e2_vec, label="FLOWERS") + axes[0].plot(orientation_vec, AEP_FLORIS_e2_vec, label="FLORIS") + axes[0].legend() + axes[1].hist((AEP_FLOWERS_e2_vec - AEP_FLORIS_e2_vec) / AEP_FLORIS_e2_vec * 100) + axes[1].set_xlabel("percent difference in FLOWERS w.r.t. FLORIS, (\\%)") + if save_plots: + fig.savefig("compare_AEP_vs_orientation.png", dpi=300, bbox_inches="tight") + plt.show() if show_plots else plt.close(fig) + +## save the output data +if save_data: + np.savez_compressed( + "rotational_consistency", + wd_vec=wd_vec, + AEP_FLOWERS_e1_vec=AEP_FLOWERS_e1_vec, + AEP_FLORIS_e1_vec=AEP_FLORIS_e1_vec, + orientation_vec=orientation_vec, + AEP_FLOWERS_e2_vec=AEP_FLOWERS_e2_vec, + AEP_FLORIS_e2_vec=AEP_FLORIS_e2_vec, + flowers_turbine=flowers_turbine, + ) diff --git a/test/flowers/system/rotation/test_rotational_consistency.py b/test/flowers/system/rotation/test_rotational_consistency.py new file mode 100644 index 00000000..c2ea28dc --- /dev/null +++ b/test/flowers/system/rotation/test_rotational_consistency.py @@ -0,0 +1,76 @@ +from pathlib import Path + +import numpy as np + +from . import test_rotational_workbench + + +class TestFLOWERSRotationalConsistency: + """ + these are tests to make sure that FLOWERS matches data we've compared + against FLORIS to make sure that the behavior of flowers matches what we + have previously calculated, particularly on a farm with no rotational + symmetries that could hide simple rotational mistakes + """ + + def setup_method(self): + path_rotational_file = Path(__file__).parent / "rotational_consistency.npz" + self.reference_pyrite_data = np.load(path_rotational_file, allow_pickle=True) + + def test_wind_direction_pyrite_match(self, subtests): + """test for a pyrite match with the orientation experiment""" + + # extract the saved off data + flowers_turbine = self.reference_pyrite_data[ + "flowers_turbine" + ].tolist() # needed to unpack + wd_vec = self.reference_pyrite_data["wd_vec"] + AEP_FLOWERS_vec_ref = self.reference_pyrite_data["AEP_FLOWERS_e1_vec"] + + # run FLOWERS using the shared toolset + AEP_FLOWERS_vec = np.zeros_like(wd_vec) + for idx, wd_val in enumerate(wd_vec): + AEP_FLOWERS_vec[idx] = test_rotational_workbench.run_FLOWERS( + flowers_turbine, wd_val=wd_val + ) + + # assert that the reference matches what we calculated + assert np.allclose(AEP_FLOWERS_vec, AEP_FLOWERS_vec_ref) + + def test_orientation_pyrite_match(self, subtests): + """test for a pyrite match with the orientation experiment""" + + # extract the saved off data + flowers_turbine = self.reference_pyrite_data[ + "flowers_turbine" + ].tolist() # needed to unpack + orientation_vec = self.reference_pyrite_data["orientation_vec"] + AEP_FLOWERS_vec_ref = self.reference_pyrite_data["AEP_FLOWERS_e2_vec"] + + # run FLOWERS using the shared toolset + AEP_FLOWERS_vec = np.zeros_like(orientation_vec) + for idx, orientation_val in enumerate(orientation_vec): + AEP_FLOWERS_vec[idx] = test_rotational_workbench.run_FLOWERS( + flowers_turbine, orientation=orientation_val + ) + + # assert that the reference matches what we calculated + assert np.allclose(AEP_FLOWERS_vec, AEP_FLOWERS_vec_ref) + + def test_FLORIS_wind_direction_reference_agreement(self, subtests): + + AEP_FLORIS_e1_vec = self.reference_pyrite_data["AEP_FLORIS_e1_vec"] + AEP_FLOWERS_e1_vec = self.reference_pyrite_data["AEP_FLOWERS_e1_vec"] + + # compute the R squared value and assert that it improves from our baseline + Rsquared = np.corrcoef(AEP_FLOWERS_e1_vec, AEP_FLORIS_e1_vec)[0, 1] ** 2 + assert Rsquared >= 0.95 * 0.935 # relax base value from reference data + + def test_FLORIS_orientation_reference_agreement(self, subtests): + + AEP_FLORIS_e2_vec = self.reference_pyrite_data["AEP_FLORIS_e2_vec"] + AEP_FLOWERS_e2_vec = self.reference_pyrite_data["AEP_FLOWERS_e2_vec"] + + # compute the R squared value and assert that it improves from our baseline + Rsquared = np.corrcoef(AEP_FLOWERS_e2_vec, AEP_FLORIS_e2_vec)[0, 1] ** 2 + assert Rsquared >= 0.95 * 0.874 # relax base value from reference data diff --git a/test/flowers/system/rotation/test_rotational_workbench.py b/test/flowers/system/rotation/test_rotational_workbench.py new file mode 100644 index 00000000..74140b7e --- /dev/null +++ b/test/flowers/system/rotation/test_rotational_workbench.py @@ -0,0 +1,79 @@ +import numpy as np +import pandas as pd + +import flowers + + +# create wraparound gaussian pulse resource for nearly single-directional wind +def f_windrose(wd_center, ws_center, wd_width=10.0, ws_width=5.0): + wd_raw = np.arange(0.0, 360.0, 5.0)[1:] + ws_raw = np.arange(0.0, 25.001, 1.0)[1:] + WS, WD = np.meshgrid(ws_raw, wd_raw) + FREQ = (np.exp(-(((WS - ws_center) / ws_width) ** 2))) * ( + np.exp(-(((WD - wd_center - 720.0) / wd_width) ** 2)) + + np.exp(-(((WD - wd_center - 360.0) / wd_width) ** 2)) + + np.exp(-(((WD - wd_center + 0.0) / wd_width) ** 2)) + + np.exp(-(((WD - wd_center + 360.0) / wd_width) ** 2)) + + np.exp(-(((WD - wd_center + 720.0) / wd_width) ** 2)) + ) + FREQ = np.maximum(FREQ, 0.001) + FREQ /= np.sum(FREQ) + df_wr = pd.DataFrame( + { + "wd": WD.flat, + "ws": WS.flat, + "freq_val": FREQ.flat, + } + ) + return (WD, WS, FREQ), df_wr + + +# generate the layout generator w/o rotational symmetries +def layout_generator( + diameter=1.0, + orientation=0.0, + factor=10 / 6, + spacing_base=5.0, + smash_over=True, +): + spacing_index_x = spacing_base * np.arange(-2, 2.1) + spacing_index_y = spacing_base * np.arange(-2, 2.1) + max_idx = np.max(np.abs(spacing_index_y)) + x_out_list = [] + y_out_list = [] + for idx in spacing_index_y: + here_factor = factor ** (idx / max_idx) + x_candidate = here_factor * spacing_index_x + if smash_over: + x_candidate = x_candidate - np.min(x_candidate) + x_candidate = x_candidate + np.min(spacing_index_x) + x_out_list.append(x_candidate) + y_out_list.append(idx * np.ones_like(spacing_index_x)) + + X_pre = diameter * np.array(x_out_list).flatten() + Y_pre = diameter * np.array(y_out_list).flatten() + + X_post, Y_post = np.array( + [ + [np.cos(np.radians(orientation)), np.sin(np.radians(orientation))], + [-np.sin(np.radians(orientation)), np.cos(np.radians(orientation))], + ] + ) @ np.vstack([X_pre, Y_pre]) + + return X_post, Y_post + + +def run_FLOWERS(flowers_turbine, wd_val=0.0, orientation=0.0): + X, Y = layout_generator(flowers_turbine["D"], orientation) + + fm = flowers.FlowersModel( + wind_rose=f_windrose( + wd_center=wd_val, + ws_center=8.0, + )[1], + layout_x=X, + layout_y=Y, + turbine=flowers_turbine, + ) + + return fm.calculate_aep() / 1e9 # return in GWh diff --git a/test/flowers/system/test_flowers_stack.py b/test/flowers/system/test_flowers_stack.py new file mode 100644 index 00000000..66500979 --- /dev/null +++ b/test/flowers/system/test_flowers_stack.py @@ -0,0 +1,109 @@ +import os +from pathlib import Path + +import numpy as np +import numpy.linalg as nLA + +import matplotlib.pyplot as plt + +import flowers +from flowers.flowers_model import FlowersModel +import flowers.tools +import flowers.visualization + +import pytest + + +def test_flowers_stack_nrel_5MW(subtests, show_plots=False): + """ + a test that puts a bunch of the pieces of FLOWERS together and makes sure + they match pyrite values + """ + + N_turb = 256 + seed_val = int.from_bytes(b"0xAR") # lock the seed for consistent results + + boundaries = [ + (-5000.0, 0.0), + (0.0, 5000.0), + (5000.0, 0.0), + (0.0, -5000.0), + ] # a big diamond grid + + # change the cwd + CWD = os.getcwd() + path_wd = Path(__file__).parent + os.chdir(path_wd) + + # load and (optionally) show the wind rose + df_wr = flowers.tools.load_wind_rose(4) + flowers.visualization.plot_wind_rose(df_wr) + plt.show() if show_plots else plt.close("all") + + # change the cwd back + os.chdir(CWD) + + # generate and (optionally) show the layout + discrete_layout = flowers.tools.discrete_layout(n_turb=N_turb, seed_val=seed_val) + flowers.visualization.plot_layout(discrete_layout[0], discrete_layout[1]) + + # run FLOWERS to get AEP + flowers_model_discrete = FlowersModel( + wind_rose=df_wr, + layout_x=discrete_layout[0], + layout_y=discrete_layout[1], + turbine="nrel_5MW", + ) + AEP_d, dAEP_d = [ + v / 1.0e9 for v in flowers_model_discrete.calculate_aep(gradient=True) + ] + AEP_d_ref = 5444.8491956346525 # pyrite value generated 22 Jan 2026 + norm_dAEP_d_ref = 0.109870151788895 # pyrite value generated 22 Jan 2026 + print(f"discrete farm AEP: {AEP_d} GW") + print(f"\tderivative magnitude: {nLA.norm(dAEP_d)} GW/m") + plt.show() if show_plots else plt.close("all") + + with subtests.test("discrete farm AEP"): + assert np.isclose(AEP_d, AEP_d_ref) + with subtests.test("discrete farm dAEP norm"): + assert np.isclose(nLA.norm(dAEP_d), norm_dAEP_d_ref) + + random_layout = flowers.tools.random_layout( + boundaries=boundaries, n_turb=N_turb, seed_val=seed_val + ) + flowers.visualization.plot_layout(random_layout[0], random_layout[1]) + + flowers_model_random = FlowersModel( + wind_rose=df_wr, + layout_x=random_layout[0], + layout_y=random_layout[1], + turbine="nrel_5MW", + ) + AEP_r, dAEP_r = [ + v / 1.0e9 for v in flowers_model_random.calculate_aep(gradient=True) + ] + AEP_r_ref = 1726.9377768748438 # pyrite value generated 22 Jan 2026 + norm_dAEP_r_ref = 0.16130660604162644 # pyrite value generated 22 Jan 2026 + print(f"random farm AEP: {AEP_r}") + print(f"\tderivative magnitude: {nLA.norm(dAEP_r)} GW/m") + plt.show() if show_plots else plt.close("all") + + with subtests.test("random farm AEP"): + assert np.isclose(AEP_r, AEP_r_ref) + with subtests.test("random farm dAEP norm"): + assert np.isclose(nLA.norm(dAEP_r), norm_dAEP_r_ref) + + +if __name__ == "__main__": + + class DummyContext: + def __enter__(self): + return self + + def __exit__(self, exc_type, exc_value, traceback): + pass + + def test(str): + return DummyContext() + + test_flowers_stack_nrel_5MW(DummyContext, show_plots=True) diff --git a/test/flowers/system/wind_roses/wr4.p b/test/flowers/system/wind_roses/wr4.p new file mode 100644 index 00000000..000a9b71 Binary files /dev/null and b/test/flowers/system/wind_roses/wr4.p differ diff --git a/test/flowers/unit/layouts/garbage5.p b/test/flowers/unit/layouts/garbage5.p new file mode 100644 index 00000000..2c72e4ec Binary files /dev/null and b/test/flowers/unit/layouts/garbage5.p differ diff --git a/test/flowers/unit/test_flowers_model.py b/test/flowers/unit/test_flowers_model.py new file mode 100644 index 00000000..763b6563 --- /dev/null +++ b/test/flowers/unit/test_flowers_model.py @@ -0,0 +1,213 @@ +from pathlib import Path + +import numpy as np +import pandas as pd + +import flowers +import ard + +import pytest + + +class TestFlowersModel: + def setup_method(self): + + self.wind_rose = pd.read_csv( + Path(flowers.__file__).parents[1] + / "examples" + / "flowers" + / "data" + / "HKW_wind_rose.csv" + ) + + self.turbine_dict = { + "D": 130.0, + # "U": None, + # "ct": None, + # "u_ct": None, + # "cp": None, + # "u_cp": None, + } + + self.layout_x, self.layout_y = [ + (self.turbine_dict["D"] * v).flat + for v in np.meshgrid( + 5 * np.arange(-3, 3.001, 1), + 5 * np.arange(-3, 3.001, 1), + ) + ] + + self.flowers_model = flowers.FlowersModel( + self.wind_rose, + layout_x=self.layout_x, + layout_y=self.layout_y, + num_terms=50, + k=0.05, + turbine="nrel_5MW", + ) + + def test_reference_NREL5MW(self, subtests): + + with subtests.test("initialization of internal values: k"): + assert np.isclose(self.flowers_model.k, 0.05) + with subtests.test("initialization of internal values: num_modes"): + assert self.flowers_model.get_num_modes() <= 50 + + # compute the AEP and compare it to a pyrite-standard value + with subtests.test("pyrite AEP value"): + AEP_ref = 904.2681563494169 # GWh + AEP_calculated = self.flowers_model.calculate_aep() / 1e9 # GWh + assert np.isclose(AEP_calculated, AEP_ref) + + def test_reinit_reference_NREL5MW(self, subtests): + + rotate_angle = np.radians(15.0) # angle to rotate layout + layout_x = ( + np.cos(rotate_angle) * self.layout_x - np.sin(rotate_angle) * self.layout_y + ) + layout_y = ( + np.sin(rotate_angle) * self.layout_x + np.cos(rotate_angle) * self.layout_y + ) + + self.flowers_model.reinitialize( + self.wind_rose, + layout_x=layout_x, + layout_y=layout_y, + num_terms=25, + k=0.075, + ) + + with subtests.test("initialization of internal values: k"): + assert np.isclose(self.flowers_model.k, 0.075) + with subtests.test("initialization of internal values: num_modes"): + assert self.flowers_model.get_num_modes() <= 25 + + # compute the AEP and compare it to a pyrite-standard value + with subtests.test("pyrite AEP value"): + AEP_ref = 976.2494933660288 # GWh + AEP_calculated = self.flowers_model.calculate_aep() / 1e9 # GWh + assert np.isclose(AEP_calculated, AEP_ref) + + def test_raise_error_reference_IEA22MW(self, subtests): + with pytest.raises(NotImplementedError): + _ = flowers.FlowersModel( + self.wind_rose, + layout_x=self.layout_x, + layout_y=self.layout_y, + num_terms=50, + k=0.05, + turbine="iea_22MW", + ) + + +class TestCustomFlowersModel: + def setup_method(self): + self.IEA_3p4_data = pd.read_csv( + Path(ard.__file__).parents[1] + / "test" + / "ard" + / "data" + / "power_thrust_table_ccblade_IEA-3p4-130-RWT.csv", + names=["u", "cp", "ct"], + ) + + self.wind_rose = pd.read_csv( + Path(flowers.__file__).parents[1] + / "examples" + / "flowers" + / "data" + / "HKW_wind_rose.csv" + ) + + self.turbine_dict = { + "D": 130.0, + "U": 25.0, + "ct": self.IEA_3p4_data["ct"].to_numpy(), + "u_ct": self.IEA_3p4_data["u"].to_numpy(), + "cp": self.IEA_3p4_data["cp"].to_numpy(), + "u_cp": self.IEA_3p4_data["u"].to_numpy(), + } + + self.layout_x, self.layout_y = [ + (self.turbine_dict["D"] * v).flat + for v in np.meshgrid( + 5 * np.arange(-3, 3.001, 1), + 5 * np.arange(-3, 3.001, 1), + ) + ] + + self.flowers_model = flowers.FlowersModel( + self.wind_rose, + layout_x=self.layout_x, + layout_y=self.layout_y, + num_terms=50, + k=0.05, + turbine=self.turbine_dict, + ) + + def test_custom_IEA3p4(self, subtests): + + # get the computed estimate + AEP_calculated = self.flowers_model.calculate_aep() / 1e9 # GWh + + # compute the AEP and compare it to a pyrite-standard value + with subtests.test("matching pyrite AEP value"): + AEP_ref = 922.0285639699603 # GWh + assert np.isclose(AEP_calculated, AEP_ref) + + # compute the capacity factor and make sure it's plausible + with subtests.test("plausible capacity factor"): + rated_production = 3.4 * 8760 * len(self.layout_x) / 1000.0 # to GWh + capacity_factor = AEP_calculated / rated_production + is_plausible_cf = (capacity_factor > 0.25) & (capacity_factor < 0.75) + assert is_plausible_cf + + def test_custom_IEA3p4_derivatives(self, subtests): + + # get the computed estimate + AEP_calculated, dAEP_calculated = [ + v / 1.0e9 for v in self.flowers_model.calculate_aep(gradient=True) + ] + + # compute the AEP through the derivative construct and pyrite-verify it + with subtests.test("matching pyrite AEP value on derivative calc"): + AEP_ref = 922.0285639699603 # GWh + assert np.isclose(AEP_calculated, AEP_ref) + + # unpack the shape and set the epsilon + Nt, Nd = dAEP_calculated.shape + eps_val = 1.0e-6 + + # break out and copy the layout variables + layout_x_orig, layout_y_orig = self.flowers_model.get_layout() + layout_x_orig = layout_x_orig.copy() + layout_y_orig = layout_y_orig.copy() + + dAEP_ping = np.zeros_like(dAEP_calculated) + # loop over turbines + for i_t in range(Nt): + for i_d in range(Nd): + # copy the original layout variables + layout_x = layout_x_orig.copy() + layout_y = layout_y_orig.copy() + # assign a ping to the right direction + layout_x[i_t] += (i_d == 0) * eps_val * np.mean(np.abs(layout_x_orig)) + layout_y[i_t] += (i_d == 1) * eps_val * np.mean(np.abs(layout_y_orig)) + self.flowers_model.reinitialize(layout_x=layout_x, layout_y=layout_y) + # compute the new AEP, the change, the derivative + AEP_plus = self.flowers_model.calculate_aep() / 1e9 + dAEP = AEP_plus - AEP_calculated + dXY = ( + (layout_y[i_t] - layout_y_orig[i_t]) + if i_d + else (layout_x[i_t] - layout_x_orig[i_t]) + ) + dAEP_ping[i_t, i_d] = dAEP / dXY + + # make sure the analytic derivative matches the ping test data + assert np.allclose( + dAEP_calculated, + dAEP_ping, + rtol=1.0e-4, + atol=1.0e-4, + ) diff --git a/test/flowers/unit/test_tool_discrete_pyrite.npz b/test/flowers/unit/test_tool_discrete_pyrite.npz new file mode 100644 index 00000000..9acd765e Binary files /dev/null and b/test/flowers/unit/test_tool_discrete_pyrite.npz differ diff --git a/test/flowers/unit/test_tool_random_pyrite.npz b/test/flowers/unit/test_tool_random_pyrite.npz new file mode 100644 index 00000000..f4f85c86 Binary files /dev/null and b/test/flowers/unit/test_tool_random_pyrite.npz differ diff --git a/test/flowers/unit/test_tools.py b/test/flowers/unit/test_tools.py new file mode 100644 index 00000000..9ce0b0a4 --- /dev/null +++ b/test/flowers/unit/test_tools.py @@ -0,0 +1,316 @@ +import os +from pathlib import Path +import pickle + +import numpy as np +import pandas as pd + +from ard.utils.test_utils import pyrite_validator +import flowers.tools + +import pytest + + +@pytest.mark.usefixtures("subtests") +def test_random_layout(subtests): + + boundaries_values = [ + (1000.0, 1000.0), + (-1000.0, 1000.0), + (-1000.0, -1000.0), + (1000.0, -1000.0), + ] # a diamond boundary + n_turb = 5 + D = 127.0 # m + min_dist = 2.5 # D_rotor + seed_val = int.from_bytes(b"0x1e") # create a big fixed integer for a random seed + + # make sure the non-specified boundary raises an error + with subtests.test("no boundary error"): + with pytest.raises(ValueError): + flowers.tools.random_layout([]) + + # make sure the non-specified turbines raises an error + with subtests.test("no turbines error"): + with pytest.raises(ValueError): + flowers.tools.random_layout(boundaries_values, n_turb=0) + + # make sure you get a boundary-conforming result w/ defaults on and no seed + xx_defaults_noseeds, yy_defaults_noseeds = flowers.tools.random_layout( + boundaries=boundaries_values, n_turb=n_turb + ) + with subtests.test("random_layout default noseed boundary conforming"): + boundary_conforming = np.all( + (xx_defaults_noseeds >= -1000.0) + & (xx_defaults_noseeds <= 1000.0) + & (yy_defaults_noseeds >= -1000.0) + & (yy_defaults_noseeds < 1000.0) + ) + assert boundary_conforming + + # make sure you get a good result w/ custom spec and no seed + xx_spec_noseeds, yy_spec_noseeds = flowers.tools.random_layout( + boundaries=boundaries_values, n_turb=n_turb, D=D, min_dist=min_dist + ) + # both boundary conforming... + with subtests.test("random_layout specified noseed boundary conforming"): + boundary_conforming = np.all( + (xx_spec_noseeds >= -1000.0) + & (xx_spec_noseeds <= 1000.0) + & (yy_spec_noseeds >= -1000.0) + & (yy_spec_noseeds < 1000.0) + ) + assert boundary_conforming + # ... and constraint satisfying + with subtests.test("random_layout specified noseed constraint satisfying"): + dist_mtx = np.sqrt( + (xx_spec_noseeds - np.atleast_2d(xx_spec_noseeds).T) ** 2 + + (yy_spec_noseeds - np.atleast_2d(yy_spec_noseeds).T) ** 2 + ) + np.fill_diagonal(dist_mtx, np.inf) + assert np.all(dist_mtx >= 2 * D) + + with subtests.test("random_layout specified seed pyrite data match"): + xx_spec_seed, yy_spec_seed = flowers.tools.random_layout( + boundaries=boundaries_values, + n_turb=n_turb, + D=D, + min_dist=min_dist, + seed_val=seed_val, + ) + dump_pickle = False + if dump_pickle: + pickle.dump( + ( + np.array(xx_spec_seed), + np.array(yy_spec_seed), + np.array(boundaries_values), + ), + open(Path(__file__).parent / "layouts" / "garbage5.p", "wb"), + ) + pyrite_validator( + data_for_validation={ + "xx": xx_spec_seed, + "yy": yy_spec_seed, + "boundaries": boundaries_values, + }, + filename_pyrite=Path(__file__).parent / "test_tool_random_pyrite", + rewrite=False, + ) + + +@pytest.mark.usefixtures("subtests") +def test_discrete_layout(subtests): + + n_turb = 5 + D = 127.0 # m + min_dist = 2.5 # D_rotor + seed_val = int.from_bytes(b"0xF5") # create a big fixed integer for a random seed + + # make sure the non-specified turbines raises an error + with subtests.test("no turbines error"): + with pytest.raises(ValueError): + flowers.tools.discrete_layout(n_turb=0) + + with subtests.test("discrete_layout specified seed pyrite data match"): + xx_spec_seed, yy_spec_seed = flowers.tools.discrete_layout( + n_turb=n_turb, + D=D, + min_dist=min_dist, + seed_val=seed_val, + ) + pyrite_validator( + data_for_validation={"xx": xx_spec_seed, "yy": yy_spec_seed}, + filename_pyrite=Path(__file__).parent / "test_tool_discrete_pyrite", + rewrite=False, + ) + + +def test_load_layout(subtests): + + # specify some figures + idx = 5 + case = "garbage" + + # change the cwd + CWD = os.getcwd() + path_wd = Path(__file__).parent + os.chdir(path_wd) + + # load the layout + layout_x, layout_y, boundaries = flowers.tools.load_layout( + idx, case, boundaries=True + ) + + # change the cwd back + os.chdir(CWD) + + # validate against pyrite data + pyrite_validator( + data_for_validation={"xx": layout_x, "yy": layout_y, "boundaries": boundaries}, + filename_pyrite=Path(__file__).parent / "test_tool_random_pyrite", + rewrite=False, + ) + + +class TestToolsWindRose: + + def setup_method(self): + self.path_csv = ( + Path(__file__).parents[3] + / "examples" + / "flowers" + / "data" + / "HKW_wind_rose.csv" + ) + self.df_wr = pd.read_csv(self.path_csv) + self.len_df_wr = len(self.df_wr) + + self.wr_idx = 4 + + def test_load_wind_rose(self, subtests): + + # change the cwd + CWD = os.getcwd() + path_wd = Path(__file__).parent + os.chdir(path_wd) + + # use the wind rose loader + df_wr_loaded = flowers.tools.load_wind_rose(self.wr_idx) + + # change the cwd back + os.chdir(CWD) + + # make sure allclose on each column + for column in self.df_wr.columns: + with subtests.test(f"loaded df matches on {column}"): + assert np.allclose(df_wr_loaded[column], self.df_wr[column]) + + def test_resample_wind_direction(self, subtests): + + # resample the wind direction + # wd_unique_orig = self.df_wr.wd.unique() + df_resampled = flowers.tools.resample_wind_direction( + self.df_wr, wd=np.arange(0, 360, 10.0) + ) + wd_unique_new = df_resampled.wd.unique() + + # get and sort the values for comparison + df_check = self.df_wr[self.df_wr.wd.isin(wd_unique_new)].sort_values( + ["ws", "wd"] + ) + df_resampled = df_resampled.sort_values(["ws", "wd"]) + + # make sure allclose on each column that's not frequency + for column in df_check: + if column == "freq_val": + continue # frequency is re-computed + with subtests.test(f"resampled df matches on {column}"): + assert np.allclose(df_resampled[column], df_check[column]) + + # resample should be approximately one + with subtests.test(f"ensure freq_val sums to one"): + assert np.isclose(np.sum(df_resampled.freq_val), 1.0, atol=0.025) + + def test_resample_wind_speed(self, subtests): + + # resample the wind speed + # wd_unique_orig = self.df_wr.wd.unique() + df_resampled = flowers.tools.resample_wind_speed( + self.df_wr, ws=np.arange(0, 25.0, 2.0) + ) + ws_unique_new = df_resampled.ws.unique() + + # get and sort the values for comparison + df_check = self.df_wr[self.df_wr.ws.isin(ws_unique_new)].sort_values( + ["ws", "wd"] + ) + df_resampled = df_resampled.sort_values(["ws", "wd"]) + + # make sure allclose on each column that's not frequency + for column in df_check: + if column == "freq_val": + continue # frequency is re-computed + with subtests.test(f"resampled df matches on {column}"): + assert np.allclose(df_resampled[column], df_check[column]) + + # resample should be approximately one + with subtests.test(f"ensure freq_val sums to one"): + assert np.isclose(np.sum(df_resampled.freq_val), 1.0, atol=0.025) + + def test_resample_average_ws_by_wd(self, subtests): + + # get the resample + df_resampled = flowers.tools.resample_average_ws_by_wd(self.df_wr) + + # check wd values + with subtests.test(f"resampling has same wd values"): + assert np.allclose(df_resampled.wd.unique(), self.df_wr.wd.unique()) + + # check freq_val + with subtests.test("frequency bins correct"): + assert np.allclose( + df_resampled.freq_val, + self.df_wr.groupby("wd").agg({"freq_val": np.sum}).freq_val, + ) + + # check average ws + with subtests.test("average speed correct"): + weighted_ws = self.df_wr.groupby("wd").apply( + lambda group: np.sum(group["ws"] * group["freq_val"]) + / np.sum(group["freq_val"]) + ) + assert np.allclose(df_resampled.ws, weighted_ws) + + +class TestToolsLookup: + + def setup_method(self): + self.u_samples = [2.0, 5.0, 8.0, 12.0, 20.0, 28.0] # check every region + self.sample_code = { + "regI": 0, + "regII_lo": 1, + "regII_hi": 2, + "regIII_lo": 3, + "regIII_hi": 4, + "regIV": 5, + } + + def test_notimplemented_raised(self, subtests): + with subtests.test("lookup cp invalid turbine raises error"): + with pytest.raises(NotImplementedError): + flowers.tools.cp_lookup( + self.u_samples, + turbine_type="the_greatest_turbine", + # cp=None, # defaulted to None + ) + with subtests.test("lookup ct invalid turbine raises error"): + with pytest.raises(NotImplementedError): + flowers.tools.ct_lookup( + self.u_samples, + turbine_type="the_greatest_turbine", + # ct=None, # defaulted to None + ) + + def test_NREL5MW(self, subtests): + + cpmax_NREL5MW = 0.436845 + ctmax_NREL5MW = 0.99 + + cp_samples = flowers.tools.cp_lookup(self.u_samples, turbine_type="nrel_5MW") + ct_samples = flowers.tools.ct_lookup(self.u_samples, turbine_type="nrel_5MW") + + with subtests.test("approx reference max cp"): + np.isclose(np.max(cp_samples), cpmax_NREL5MW) + with subtests.test("approx reference max ct"): + np.isclose(np.max(ct_samples), ctmax_NREL5MW) + + with subtests.test("region I cp approx zero"): + np.isclose(cp_samples[self.sample_code["regI"]], 0.0) + with subtests.test("region I ct approx zero"): + np.isclose(ct_samples[self.sample_code["regI"]], 0.0) + + with subtests.test("region IV cp approx zero"): + np.isclose(cp_samples[self.sample_code["regIV"]], 0.0) + with subtests.test("region IV ct approx zero"): + np.isclose(ct_samples[self.sample_code["regIV"]], 0.0) diff --git a/test/flowers/unit/wind_roses/wr4.p b/test/flowers/unit/wind_roses/wr4.p new file mode 100644 index 00000000..000a9b71 Binary files /dev/null and b/test/flowers/unit/wind_roses/wr4.p differ diff --git a/test/run_local_test_system.sh b/test/run_local_test_system.sh index 28729015..a2c50f8d 100644 --- a/test/run_local_test_system.sh +++ b/test/run_local_test_system.sh @@ -1,10 +1,32 @@ #!/bin/bash -if python -c "import optiwindnet" 2>/dev/null ; then - pytest --cov=ard --cov-report=html test/ard/system -else - pytest --cov=ard --cov-report=html test/ard/system --cov-config=.coveragerc_no_optiwindnet -fi -rm -rf test/ard/system/layout/problem*_out +# use case to branch and run cases +case "$1" in + # first branch: one subpackage at a time: valid choices separated by pipe + ard|flowers) + CASES_TO_RUN=$1 # select-your-own + ;; + # default (no arg) behavior: run all + ""|all) + CASES_TO_RUN="ard flowers" # default case ("all") -> ard then flowers + ;; + # error case: print error message and return >0 + *) + echo "$1 NOT FOUND IN VALID CASES!" # error case + return 999 # + ;; +esac -# +# cycle over indexed cases +IDX_CASE=0 +for CASE in $CASES_TO_RUN; do + FLAGS="--cov-report=html" # default flag: html coverage report + if [[ $IDX_CASE -gt 0 ]] ; then + FLAGS="$FLAGS --cov-append" # append after the first case for multiple + fi + if [[ "$CASE" == "ard" ]] && ! python -c "import optiwindnet" 2>/dev/null ; then + FLAGS="$FLAGS --cov-config=.coveragerc_no_optiwindnet" # special case + fi + pytest --cov=$CASE ${FLAGS} test/$CASE/system # run the case + ((IDX_CASE++)) +done diff --git a/test/run_local_test_unit.sh b/test/run_local_test_unit.sh index a64346a7..c8b7553d 100644 --- a/test/run_local_test_unit.sh +++ b/test/run_local_test_unit.sh @@ -1,10 +1,32 @@ #!/bin/bash -if python -c "import optiwindnet" 2>/dev/null ; then - pytest --cov=ard --cov-report=html test/ard/unit -else - pytest --cov=ard --cov-report=html test/ard/unit --cov-config=.coveragerc_no_optiwindnet -fi -rm -rf test/ard/unit/layout/problem*_out +# use case to branch and run cases +case "$1" in + # first branch: one subpackage at a time: valid choices separated by pipe + ard|flowers) + CASES_TO_RUN=$1 # select-your-own + ;; + # default (no arg) behavior: run all + ""|all) + CASES_TO_RUN="ard flowers" # default case ("all") -> ard then flowers + ;; + # error case: print error message and return >0 + *) + echo "$1 NOT FOUND IN VALID CASES!" # error case + return 999 # + ;; +esac -# +# cycle over indexed cases +IDX_CASE=0 +for CASE in $CASES_TO_RUN; do + FLAGS="--cov-report=html" # default flag: html coverage report + if [[ $IDX_CASE -gt 0 ]] ; then + FLAGS="$FLAGS --cov-append" # append after the first case for multiple + fi + if [[ "$CASE" == "ard" ]] && ! python -c "import optiwindnet" 2>/dev/null ; then + FLAGS="$FLAGS --cov-config=.coveragerc_no_optiwindnet" # special case + fi + pytest --cov=$CASE ${FLAGS} test/$CASE/unit # run the case + ((IDX_CASE++)) +done