diff --git a/.github/workflows/continuous-integration-workflow.yaml b/.github/workflows/continuous-integration-workflow.yaml index 8159a0ba..88a187eb 100644 --- a/.github/workflows/continuous-integration-workflow.yaml +++ b/.github/workflows/continuous-integration-workflow.yaml @@ -24,6 +24,7 @@ jobs: pip install -e ".[develop]" pip install git+https://github.com/NREL/electrolyzer.git pip install https://github.com/NREL/SEAS/blob/main/SEAS.tar.gz?raw=true + pip install git+https://github.com/NREL/floris.git@v4 # - uses: pre-commit/action@v3.0.0 - name: Run ruff run: | diff --git a/.gitignore b/.gitignore index f3bcfa89..2c56ddae 100644 --- a/.gitignore +++ b/.gitignore @@ -34,12 +34,14 @@ hercules/local_amr_wind_demo/sample_copy.nc #Ignore csv files *.csv +!tests/test_inputs/amr_standin_data.csv # Some output files to ignore t_00* logdummy loghercules logstandin +logfloris *echo *out-example.json diff --git a/docs/install_old.md b/docs/install_old.md index 217aba97..a7015ef1 100644 --- a/docs/install_old.md +++ b/docs/install_old.md @@ -70,6 +70,13 @@ NREL's PySAM software is also required for hercules. To install, use pip install nrel-pysam==4.2.0 ``` +You may also want to run the FLORIS standin for AMR-Wind for a steady-state representation +of wind farm flows. In this case, run +``` +pip install git+https://github.com/NREL/floris.git@v4 +``` +to get the correct branch of FLORIS installed. + If you run hercules and get an error that `pyyaml` is missing, you may also need to install it using ``` conda install -c conda-forge pyyaml diff --git a/example_case_folders/08_floris_only/amr_input.inp b/example_case_folders/08_floris_only/amr_input.inp new file mode 100644 index 00000000..4b811607 --- /dev/null +++ b/example_case_folders/08_floris_only/amr_input.inp @@ -0,0 +1,165 @@ +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# SIMULATION STOP # +#.......................................# +time.stop_time = 100.0 # Max (simulated) time to evolve +time.max_step = -1 # Max number of time steps + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# TIME STEP COMPUTATION # +#.......................................# +time.fixed_dt = 0.5 # Use this constant dt if > 0 +time.cfl = 0.95 # CFL factor + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# INPUT AND OUTPUT # +#.......................................# +time.plot_interval = 3600 # Steps between plot files +time.checkpoint_interval = 3600 # Steps between checkpoint files +io.restart_file = "/projects/ssc/amr_precursors/b_abl_neutral_lowTI_redo/chk14400" + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# PHYSICS # +#.......................................# +incflo.gravity = 0. 0. -9.81 # Gravitational force (3D) +incflo.density = 1.0 # Reference density + +incflo.use_godunov = 1 +incflo.godunov_type = weno_z +incflo.diffusion_type = 1 +transport.viscosity = 1.0e-5 +transport.laminar_prandtl = 0.7 +transport.turbulent_prandtl = 0.3333 +turbulence.model = OneEqKsgsM84 + +incflo.physics = ABL Actuator +ICNS.source_terms = BoussinesqBuoyancy CoriolisForcing ABLMeanBoussinesq ActuatorForcing +TKE.source_terms = KsgsM84Src +BoussinesqBuoyancy.reference_temperature = 300.0 +CoriolisForcing.latitude = 41.3 +ABLForcing.abl_forcing_height = 90 +incflo.velocity = 6.928203230275509 4.0 0.0 + + +# Atmospheric boundary layer +ABL.temperature_heights = 0.0 700.0 800.0 1280.0 +ABL.temperature_values = 300.0 300.0 308.0 309.44 +ABL.reference_temperature = 300.0 +ABL.kappa = .40 +ABL.surface_roughness_z0 = 1.0E-4 +ABL.Uperiods = 25.0 +ABL.Vperiods = 25.0 +ABL.cutoff_height = 50.0 +ABL.deltaU = 1.0 +ABL.deltaV = 1.0 +ABL.normal_direction = 2 +ABL.perturb_ref_height = 50.0 +ABL.perturb_temperature = false +ABL.perturb_velocity = true +ABL.stats_output_format = netcdf +ABL.stats_output_frequency = 1 +ABL.surface_temp_flux = 0.00 +ABL.wall_shear_stress_type = "Moeng" + +ABL.bndry_file = "/projects/ssc/amr_precursors/b_abl_neutral_lowTI_redo/bndry_files" +ABL.bndry_io_mode = 1 +ABL.bndry_planes = ylo xlo # I'm (Paul) adding this but not sure if I have to +ABL.bndry_var_names = velocity temperature tke + + +# Output boundary files +ABL.bndry_planes = ylo xlo +ABL.bndry_output_start_time = 7200.0 +ABL.bndry_var_names = velocity temperature tke +ABL.bndry_output_format = native +ABL.stats_output_frequency = 1 +ABL.stats_output_format = netcdf + +# Whether to use helics +helics.activated = true +helics.broker_port = 32000 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# ADAPTIVE MESH REFINEMENT # +#.......................................# +amr.n_cell = 512 512 128 # Grid cells at coarsest AMRlevel +amr.max_level = 0 # Max AMR level in hierarchy + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# GEOMETRY # +#.......................................# +geometry.prob_lo = 0. 0. 0. # Lo corner coordinates +geometry.prob_hi = 5120. 5120. 1280. # Hi corner coordinates +geometry.is_periodic = 0 0 0 +xlo.type = "mass_inflow" +xlo.density = 1.0 +xlo.temperature = 0.0 # value required but ignored +xlo.tke = 0.0 +xhi.type = "pressure_outflow" + +ylo.type = "mass_inflow" +ylo.density = 1.0 +ylo.temperature = 0.0 +ylo.tke = 0.0 +yhi.type = "pressure_outflow" + +# Boundary conditions +zlo.type = "wall_model" +zlo.tke_type = "zero_gradient" + +zhi.type = "slip_wall" +zhi.temperature_type = "fixed_gradient" +zhi.temperature = 0.003 # tracer is used to specify potential temperature gradient + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# VERBOSITY # +#.......................................# +incflo.verbose = 0 # incflo_level + + + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# SAMPLING # +#.......................................# +incflo.post_processing = samplingPlane samplingLine + +samplingPlane.output_frequency = 600 +samplingPlane.labels = z_plane +samplingPlane.fields = velocity temperature +samplingPlane.z_plane.type = PlaneSampler +samplingPlane.z_plane.axis1 = 5110 0.0 0.0 +samplingPlane.z_plane.axis2 = 0.0 5110 0.0 +samplingPlane.z_plane.origin = 5.0 5.0 0.0 +samplingPlane.z_plane.num_points = 512 512 +samplingPlane.z_plane.normal = 0.0 0.0 1.0 +samplingPlane.z_plane.offsets = 5.0 85.0 155.0 255.0 + + +samplingLine.output_frequency = 1 +samplingLine.labels = z_line +samplingLine.fields = velocity temperature +samplingLine.z_line.type = LineSampler +samplingLine.z_line.num_points = 128 +samplingLine.z_line.start = 5.0 5.0 5.0 +samplingLine.z_line.end = 5.0 5.0 1275.0 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# TURBINES # +#.......................................# + +Actuator.type = UniformCtDisk +Actuator.UniformCtDisk.rotor_diameter = 126.0 +Actuator.UniformCtDisk.hub_height = 90.0 +Actuator.UniformCtDisk.thrust_coeff = 0.0 0.0 1.132034888 0.999470963 0.917697381 0.860849503 0.815371198 0.811614904 0.807939328 0.80443352 0.800993851 0.79768116 0.794529244 0.791495834 0.788560434 0.787217182 0.787127977 0.785839257 0.783812219 0.783568108 0.783328285 0.781194418 0.777292539 0.773464375 0.769690236 0.766001924 0.762348072 0.758760824 0.755242872 0.751792927 0.748434131 0.745113997 0.717806682 0.672204789 0.63831272 0.610176496 0.585456847 0.563222111 0.542912273 0.399312061 0.310517829 0.248633226 0.203543725 0.169616419 0.143478955 0.122938861 0.106515296 0.093026095 0.081648606 0.072197368 0.064388275 0.057782745 0.0 0.0 +Actuator.UniformCtDisk.wind_speed = 0.0 2.9 3.0 4.0 5.0 6.0 7.0 7.1 7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9 8.0 9.0 10.0 10.1 10.2 10.3 10.4 10.5 10.6 10.7 10.8 10.9 11.0 11.1 11.2 11.3 11.4 11.5 11.6 11.7 11.8 11.9 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 25.1 50.0 +Actuator.UniformCtDisk.epsilon = 10.0 +Actuator.UniformCtDisk.density = 1.225 +Actuator.UniformCtDisk.diameters_to_sample = 1.0 +Actuator.UniformCtDisk.num_points_r = 20 +Actuator.UniformCtDisk.num_points_t = 5 + +Actuator.UniformCtDisk.yaw = 240.0 + +Actuator.labels = T00 T01 + +Actuator.T00.base_position = 2000.0 2000.0 0.0 +Actuator.T01.base_position = 2500.0 2500.0 0.0 diff --git a/example_case_folders/08_floris_only/bash_script.sh b/example_case_folders/08_floris_only/bash_script.sh new file mode 100755 index 00000000..9699dfbb --- /dev/null +++ b/example_case_folders/08_floris_only/bash_script.sh @@ -0,0 +1,28 @@ +# Example bash for running things locally +# I just run these one at a t time + +# A lot of modules and conda stuff +conda activate hercules + +# Set the helics port to use: +export HELICS_PORT=32000 + +#make sure you use the same port number in the amr_input.inp and hercules_input_000.yaml files. + +# Clear old log files for clarity +rm loghercules logfloris + +# Set up the helics broker +helics_broker -t zmq -f 2 --loglevel="debug" --local_port=$HELICS_PORT & +#helics_broker -f 2 --consoleloglevel=trace --loglevel=debug --local_port=$HELICS_PORT >> loghelics & + +# Need to set this to your hercules folder +# cd /home/pfleming/hercules/hercules +python3 hercules_runscript.py hercules_input_000.yaml >> loghercules 2>&1 & # Start the controller center and pass in input file + + +python3 floris_runscript.py amr_input.inp >> logfloris 2>&1 +# Now go back to scratch folder and launch the job + +# cd /scratch/pfleming/c2c/example_sim_02 +# mpirun -n 72 /home/pfleming/amr-wind/build/amr_wind amr_input.inp >> logamr diff --git a/example_case_folders/08_floris_only/floris_runscript.py b/example_case_folders/08_floris_only/floris_runscript.py new file mode 100644 index 00000000..7fcdae89 --- /dev/null +++ b/example_case_folders/08_floris_only/floris_runscript.py @@ -0,0 +1,15 @@ +import sys + +from hercules.floris_standin import launch_floris + +# Check that one command line argument was given +if len(sys.argv) != 2: + raise Exception("Usage: python floris_runscript.py ") + +# # Get the first command line argument +# This is the name of the file to read +amr_input_file = sys.argv[1] +print(f"Running FLORIS standin with input file: {amr_input_file}") + + +launch_floris(amr_input_file) diff --git a/example_case_folders/08_floris_only/hercules_input_000.yaml b/example_case_folders/08_floris_only/hercules_input_000.yaml new file mode 100644 index 00000000..592854f6 --- /dev/null +++ b/example_case_folders/08_floris_only/hercules_input_000.yaml @@ -0,0 +1,49 @@ +# Input YAML for emy_python + +# Name +name: example_000 + +### +# Describe this emulator setup +description: AMR-wind (Standin) Only + +dt: 0.5 + +hercules_comms: + + amr_wind: + + wind_farm_0: + type: amr_wind_local #options are amr_wind or amr_wind_local + amr_wind_input_file: amr_input.inp + + helics: + + config: + name: hercules # What is the purpose of this name + use_dash_frontend: False + KAFKA: False + KAFKA_topics: EMUV1py + helics: + # deltat: 1 # This will be assigned in software + subscription_topics: [status] + publication_topics: [control] + endpoints: [] + helicsport : 32000 + publication_interval: 1 + endpoint_interval: 1 + starttime: 0 + stoptime: 100 + + Agent: ControlCenter + +py_sims: + +controller: + + + + + + + diff --git a/example_case_folders/08_floris_only/hercules_runscript.py b/example_case_folders/08_floris_only/hercules_runscript.py new file mode 100644 index 00000000..a98a6073 --- /dev/null +++ b/example_case_folders/08_floris_only/hercules_runscript.py @@ -0,0 +1,17 @@ +import sys + +from hercules.controller_standin import ControllerStandin +from hercules.emulator import Emulator +from hercules.py_sims import PySims +from hercules.utilities import load_yaml + +input_dict = load_yaml(sys.argv[1]) + + +controller = ControllerStandin(input_dict) +py_sims = PySims(input_dict) + + +emulator = Emulator(controller, py_sims, input_dict) +emulator.run_helics_setup() +emulator.enter_execution(function_targets=[], function_arguments=[[]]) diff --git a/example_case_folders/08_floris_only/readme.txt b/example_case_folders/08_floris_only/readme.txt new file mode 100644 index 00000000..1a3c77fd --- /dev/null +++ b/example_case_folders/08_floris_only/readme.txt @@ -0,0 +1 @@ +This example runs FLORIS instead of AMR-Wind \ No newline at end of file diff --git a/hercules/amr_wind_standin.py b/hercules/amr_wind_standin.py index 30501fb3..439a828d 100644 --- a/hercules/amr_wind_standin.py +++ b/hercules/amr_wind_standin.py @@ -184,12 +184,16 @@ def run(self): logger.info("Calculating simulation time: %.1f" % sim_time_s) # Compute the turbine power using a simple formula + if self.message_from_server is not None: + yaw_angles = self.message_from_server[-self.num_turbines :] + else: + yaw_angles = None ( amr_wind_speed, amr_wind_direction, turbine_powers, turbine_wind_directions, - ) = self.get_step(sim_time_s) + ) = self.get_step(sim_time_s, yaw_angles) # ================================================================ # Communicate with control center @@ -232,7 +236,7 @@ def run(self): # TODO cleanup code to move publish and subscribe here. - def get_step(self, sim_time_s): + def get_step(self, sim_time_s, yaw_angles=None): """Retreive or calculate wind speed, direction, and turbine powers Input: @@ -327,7 +331,7 @@ def launch_amr_wind_standin(amr_input_file, amr_standin_data_file=None): "endpoint_interval": 1, "starttime": 0, "stoptime": temp["stop_time"], - "Agent": "dummy_amr_wind", + "Agent": "amr_wind_standin", } if amr_standin_data_file is not None: diff --git a/hercules/emulator.py b/hercules/emulator.py index e13a3a3c..b65ef23b 100644 --- a/hercules/emulator.py +++ b/hercules/emulator.py @@ -371,6 +371,9 @@ def read_amr_wind_input(self, amr_wind_input): if "Actuator.labels" in line: turbine_labels = line.split()[2:] num_turbines = len(turbine_labels) + for line in Lines: + if "Actuator.type" in line: + actuator_type = line.split()[-1] self.num_turbines = num_turbines print("Number of turbines in amrwind: ", num_turbines) @@ -397,7 +400,7 @@ def read_amr_wind_input(self, amr_wind_input): # Find the diameter for line in Lines: - if "rotor_diameter" in line: + if "Actuator.%s.rotor_diameter" % actuator_type in line: D = float(line.split()[-1]) # Get the turbine locations diff --git a/hercules/floris_standin.py b/hercules/floris_standin.py new file mode 100644 index 00000000..288816fe --- /dev/null +++ b/hercules/floris_standin.py @@ -0,0 +1,304 @@ +# Copyright 2022 NREL + +# Licensed under the Apache License, Version 2.0 (the "License"); you may not +# use this file except in compliance with the License. You may obtain a copy of +# the License at http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, WITHOUT +# WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the +# License for the specific language governing permissions and limitations under +# the License. + +# This script implements a test client to test out the server against +# It is based on code from +# https://github.com/TUDelft-DataDrivenControl/SOWFA/blob/master/exampleCases/example.12.piso.NREL5MW.ADM.zmqSSC.python/ssc/testclient.py + +# The basic operation for this script is to pretend to be the wind farm simulator by: +# - First connecting to the front-end server +# - Than in a loop: +# - - Send the measurement values of 4 turbines +# - - Receive the wind speed and wind direction measurements +# - - Update the turbine measurements +# - - Sleep for 1 s + +import logging +import sys + +import numpy as np +from floris.tools import FlorisInterface +from floris.turbine_library import build_cosine_loss_turbine_dict + +from hercules.amr_wind_standin import AMRWindStandin, read_amr_wind_input + +# Set up the logger +# Useful for when running on eagle +logging.basicConfig( + level=logging.DEBUG, + format="%(asctime)s %(name)-12s %(levelname)-8s %(message)s", + datefmt="%Y-%m-%d %H:%M", + filename="log_test_client.log", + filemode="w", +) +logger = logging.getLogger("amr_wind_standin") + +# Perhaps a small hack to also send log to the terminal outout +logging.getLogger().addHandler(logging.StreamHandler(sys.stdout)) + +# Make an announcement +logger.info("Emulator amr_wind_standin (standing in for AMR-Wind) connecting to server") + + +# Define a function to read the amrwind input file +# Note simply copied from emulator +def construct_floris_from_amr_input(amr_wind_input): + # Probably want a file not found error instead + + with open(amr_wind_input) as fp: + Lines = fp.readlines() + + # Get the turbine locations + layout_x = [] + layout_y = [] + for line in Lines: + if ".base_position" in line: + loc = [float(d) for d in line.split()[2:]] + layout_x.append(loc[0]) + layout_y.append(loc[1]) + + # Check that uniform disk is specified + for line in Lines: + if "Actuator.type" in line: + acuator_type = line.split()[2] + if acuator_type == "UniformCtDisk": + pass + else: + raise NotImplementedError("FLORIS standin requires UniformCtDisk actuators.") + + # Get the turine parameters + for line in Lines: + if acuator_type + ".rotor_diameter" in line: + rotor_diameter = float(line.split()[2]) + for line in Lines: + if acuator_type + ".hub_height" in line: + hub_height = float(line.split()[2]) + for line in Lines: + if acuator_type + ".density" in line: + ref_air_density = float(line.split()[2]) + + # construct turbine thrust and power curves + for line in Lines: + if acuator_type + ".thrust_coeff" in line: + thrust_coefficient = [float(d) for d in line.split()[2:]] + for line in Lines: + if acuator_type + ".wind_speed" in line: + wind_speed = [float(d) for d in line.split()[2:]] + # The power curve needs to be constructed from available data + thrust_coefficient = np.array(thrust_coefficient) + if (thrust_coefficient < 0).any() or (thrust_coefficient > 1).any(): + print("Clipping thrust_coefficient to (0, 1) interval.") + thrust_coefficient = np.clip(thrust_coefficient, 0.0, 1.0) + ai = (1 - np.sqrt(1 - np.array(thrust_coefficient))) / 2 + power_coefficient = 4 * ai * (1 - ai) ** 2 + turbine_data_dict = { + "wind_speed": wind_speed, + "thrust_coefficient": thrust_coefficient, + "power_coefficient": list(power_coefficient), + } + + # Build the turbine dictionary as expected by FLORIS + turb_dict = build_cosine_loss_turbine_dict( + turbine_data_dict=turbine_data_dict, + turbine_name="FLORIS-standin-turbine", + hub_height=hub_height, + rotor_diameter=rotor_diameter, + ref_air_density=ref_air_density, + ) + + # load a default model + fi = FlorisInterface(default_floris_dict) + fi.reinitialize( + layout_x=layout_x, layout_y=layout_y, turbine_type=[turb_dict] * len(layout_x) + ) + + return fi + + +class FlorisStandin(AMRWindStandin): + def __init__(self, config_dict, amr_input_file, amr_standin_data_file=None): + super(FlorisStandin, self).__init__( + config_dict=config_dict, + amr_wind_input=amr_input_file, + amr_standin_data_file=amr_standin_data_file, + ) + + # Construct the floris object + self.fi = construct_floris_from_amr_input(amr_input_file) + + # Get the number of turbines + self.num_turbines = len(self.fi.layout_x) + + # Print the number of turbines + logger.info("Number of turbines: {}".format(self.num_turbines)) + + # Initialize storage + self.yaw_angles_stored = [0.0] * self.num_turbines + + def get_step(self, sim_time_s, yaw_angles=None): + """Retreive or calculate wind speed, direction, and turbine powers + + Input: + sim_time_s: simulation time step + + Output: + amr_wind_speed: wind speed at current time step + amr_wind_direction: wind direction at current time step + turbine_powers: turbine powers at current time step + """ + + if hasattr(self, "standin_data"): + amr_wind_speed = np.interp( + sim_time_s, + self.standin_data["time"], + self.standin_data["amr_wind_speed"], + ) + amr_wind_direction = np.interp( + sim_time_s, + self.standin_data["time"], + self.standin_data["amr_wind_direction"], + ) + + else: + amr_wind_speed = 8.0 + amr_wind_direction = 240.0 + + turbine_wind_directions = [amr_wind_direction] * self.num_turbines + + # Compute the turbine power using FLORIS + self.fi.reinitialize(wind_speeds=[amr_wind_speed], wind_directions=[amr_wind_direction]) + + if yaw_angles is not None: + yaw_misalignments = (amr_wind_direction - np.array(yaw_angles))[ + None, : + ] # TODO: remove 2 + + if (np.abs(yaw_misalignments) > 45).any(): # Bad yaw angles + print( + ( + "Large yaw misalignment detected. " + + "Wind direction: {0:.2f} deg, ".format(amr_wind_direction) + + "Yaw angles: " + ), + yaw_angles, + "Using previous yaw angles.", + ) + yaw_misalignments = (amr_wind_direction - np.array(self.yaw_angles_stored))[None, :] + else: # Reasonable yaw angles, save in case bad angles received later + self.yaw_angles_stored = yaw_angles + else: + yaw_misalignments = yaw_angles + self.fi.calculate_wake(yaw_angles=yaw_misalignments) + # This converts the output power from Floris (in Watts) to kW (standard for Hercules) + turbine_powers = (self.fi.get_turbine_powers() / 1000).flatten().tolist() # in kW + + return ( + amr_wind_speed, + amr_wind_direction, + turbine_powers, + turbine_wind_directions, + ) + + def process_endpoint_event(self, msg): + pass + + def process_periodic_endpoint(self): + pass + + def process_periodic_publication(self): + # Periodically publish data to the surrpogate + pass + + def process_subscription_messages(self, msg): + pass + + +def launch_floris(amr_input_file, amr_standin_data_file=None): + temp = read_amr_wind_input(amr_input_file) + + config = { + "name": "floris_standin", + "gridpack": {}, + "helics": { + "deltat": temp["dt"], + "subscription_topics": ["control"], + "publication_topics": ["status"], + "endpoints": [], + "helicsport": temp["helics_port"], + }, + "publication_interval": 1, + "endpoint_interval": 1, + "starttime": 0, + "stoptime": temp["stop_time"], + "Agent": "floris_standin", + } + + obj = FlorisStandin(config, amr_input_file, amr_standin_data_file) + + obj.run_helics_setup() + obj.enter_execution(function_targets=[], function_arguments=[[]]) + + +default_floris_dict = { + "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_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}, + }, + "wake_turbulence_parameters": { + "crespo_hernandez": {"initial": 0.1, "constant": 0.5, "ai": 0.8, "downstream": -0.32} + }, + "wake_velocity_parameters": { + "gauss": {"alpha": 0.58, "beta": 0.077, "ka": 0.38, "kb": 0.004}, + "jensen": {"we": 0.05}, + }, + }, + "farm": { + "layout_x": [0.0], + "layout_y": [0.0], + "turbine_type": ["nrel_5MW"], + }, + "flow_field": { + "wind_speeds": [8.0], + "wind_directions": [270.0], + "wind_veer": 0.0, + "wind_shear": 0.12, + "air_density": 1.225, + "turbulence_intensity": 0.06, + "reference_wind_height": 90.0, + }, + "name": "GCH_for_FlorisStandin", + "description": "FLORIS Gauss Curl Hybrid model standing in for AMR-Wind", + "floris_version": "v4.x", +} diff --git a/tests/floris_standin_test.py b/tests/floris_standin_test.py new file mode 100644 index 00000000..45c3767f --- /dev/null +++ b/tests/floris_standin_test.py @@ -0,0 +1,156 @@ +from pathlib import Path + +import numpy as np +from floris.tools import FlorisInterface +from hercules.amr_wind_standin import AMRWindStandin +from hercules.floris_standin import ( + construct_floris_from_amr_input, + default_floris_dict, + FlorisStandin, +) +from SEAS.federate_agent import FederateAgent + +AMR_INPUT = Path(__file__).resolve().parent / "test_inputs" / "amr_input_florisstandin.inp" +AMR_EXTERNAL_DATA = Path(__file__).resolve().parent / "test_inputs" / "amr_standin_data.csv" + +CONFIG = { + "name": "floris_standin", + "gridpack": {}, + "helics": { + "deltat": 1.0, + "subscription_topics": ["control"], + "publication_topics": ["status"], + "endpoints": [], + "helicsport": None, + }, + "publication_interval": 1, + "endpoint_interval": 1, + "starttime": 0, + "stoptime": 10.0, + "Agent": "floris_standin", +} + + +def test_construct_floris_from_amr_input(): + fi_test = construct_floris_from_amr_input(AMR_INPUT) + assert isinstance(fi_test, FlorisInterface) + + +def test_FlorisStandin_instantiation(): + # Check instantiates correctly + floris_standin = FlorisStandin(CONFIG, AMR_INPUT) + + # Check inheritance + assert isinstance(floris_standin, FederateAgent) + assert isinstance(floris_standin, AMRWindStandin) + + # Get FLORIS equivalent, match layout and turbines + fi_true = FlorisInterface(default_floris_dict) + fi_true.reinitialize( + layout_x=floris_standin.fi.layout_x, + layout_y=floris_standin.fi.layout_y, + turbine_type=floris_standin.fi.floris.farm.turbine_definitions, + ) + + assert fi_true.floris.as_dict() == floris_standin.fi.floris.as_dict() + + +def test_FlorisStandin_get_step(): + floris_standin = FlorisStandin(CONFIG, AMR_INPUT) + + # Get FLORIS equivalent, match layout and turbines + fi_true = FlorisInterface(default_floris_dict) + fi_true.reinitialize( + layout_x=floris_standin.fi.layout_x, + layout_y=floris_standin.fi.layout_y, + turbine_type=floris_standin.fi.floris.farm.turbine_definitions, + ) + + default_wind_direction = 240.0 # Matches default in FlorisStandin + default_wind_speed = 8.0 # Matches default in FlorisStandin + + # Test with None yaw angles + fs_ws, fs_wd, fs_tp, fs_twd = floris_standin.get_step(5.0) + fi_true.reinitialize(wind_speeds=[default_wind_speed], wind_directions=[default_wind_direction]) + fi_true.calculate_wake() + fi_true_tp = fi_true.get_turbine_powers() / 1000 # kW expected + + assert fs_ws == default_wind_speed + assert fs_wd == default_wind_direction + assert fs_twd == [default_wind_direction] * 2 + assert np.allclose(fs_tp, fi_true_tp.flatten().tolist()) + + # Test with aligned turbines + yaw_angles = [240.0, 240.0] + fs_ws, fs_wd, fs_tp, fs_twd = floris_standin.get_step(5.0, yaw_angles) + fi_true.reinitialize(wind_speeds=[default_wind_speed], wind_directions=[default_wind_direction]) + fi_true.calculate_wake() # Aligned in any case + fi_true_tp = fi_true.get_turbine_powers() / 1000 # kW expected + + assert np.allclose(fs_tp, fi_true_tp.flatten().tolist()) + + # Test with misaligned turbines + yaw_angles = [260.0, 230.0] + fs_ws, fs_wd, fs_tp, fs_twd = floris_standin.get_step(5.0, yaw_angles) + fi_true.reinitialize(wind_speeds=[default_wind_speed], wind_directions=[default_wind_direction]) + fi_true.calculate_wake() # Don't expect to work + fi_true_tp = fi_true.get_turbine_powers() + assert not np.allclose(fs_tp, fi_true_tp.flatten().tolist()) + + # Correct yaw angles + fi_true.calculate_wake(yaw_angles=default_wind_direction - np.array([yaw_angles])) + fi_true_tp = fi_true.get_turbine_powers() / 1000 # kW expected + assert np.allclose(fs_tp, fi_true_tp.flatten().tolist()) + + +def test_FlorisStandin_with_standin_data(): + floris_standin = FlorisStandin(CONFIG, AMR_INPUT, AMR_EXTERNAL_DATA) + + yaw_angles_all = [ + [240.0, 240.0], + [240.0, 240.0], # Step up from 8 to 10 m/s + [240.0, 240.0], + [240.0, 240.0], # Step back down to 8 m/s + [240.0, 240.0], # wd changes to 270.0 here + [270.0, 270.0], # change to match wd + [270.0, 270.0], + [250.0, 270.0], # Apply simple make steering + [250.0, 270.0], + [250.0, 270.0], + ] + + # Initialize storage + fs_ws_all = [] + fs_wd_all = [] + fs_tp_all = [] + fs_twd_all = [] + for i, s in enumerate(np.arange(0, 10.0, 1.0)): + fs_ws, fs_wd, fs_tp, fs_twd = floris_standin.get_step(s, yaw_angles=yaw_angles_all[i]) + fs_ws_all.append(fs_ws) + fs_wd_all.append(fs_wd) + fs_tp_all.append(fs_tp) + fs_twd_all.append(fs_twd) + + # Check standin data mapped over correctly + assert fs_ws_all == floris_standin.standin_data.amr_wind_speed.to_list() + assert fs_wd_all == floris_standin.standin_data.amr_wind_direction.to_list() + assert np.allclose( + np.array(fs_twd_all)[:, 0], floris_standin.standin_data.amr_wind_direction.values + ) + + # Check power behaves as expected + # Same condition for each + fs_tp_all = np.array(fs_tp_all) + assert np.allclose(fs_tp_all[0, :], fs_tp_all[3, :]) + # Higher power at 10 than 8 m/s + assert (np.array(fs_tp_all[0, :]) < np.array(fs_tp_all[1, :])).all() + # Same power at upstream turbine, lower power at downstream turbine when aligned + assert fs_tp_all[5, 0] == fs_tp_all[0, 0] + assert fs_tp_all[5, 1] < fs_tp_all[0, 1] + # Lower power at upstream turbine, higher power at downstream turbine when steering, + # total power uplift + assert fs_tp_all[7, 0] < fs_tp_all[6, 0] + assert fs_tp_all[7, 1] > fs_tp_all[6, 1] + assert fs_tp_all[7, :].sum() > fs_tp_all[6, :].sum() + # More power steering at 10m/s than 8m/s + assert fs_tp_all[9, :].sum() > fs_tp_all[7, :].sum() diff --git a/tests/test_inputs/amr_input_florisstandin.inp b/tests/test_inputs/amr_input_florisstandin.inp new file mode 100644 index 00000000..3de6ac64 --- /dev/null +++ b/tests/test_inputs/amr_input_florisstandin.inp @@ -0,0 +1,181 @@ +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# SIMULATION STOP # +#.......................................# +time.stop_time = 100.0 # Max (simulated) time to evolve +time.max_step = -1 # Max number of time steps + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# TIME STEP COMPUTATION # +#.......................................# +time.fixed_dt = 0.5 # Use this constant dt if > 0 +time.cfl = 0.95 # CFL factor + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# INPUT AND OUTPUT # +#.......................................# +time.plot_interval = 3600 # Steps between plot files +time.checkpoint_interval = 3600 # Steps between checkpoint files +io.restart_file = "/projects/ssc/amr_precursors/b_abl_neutral_lowTI_redo/chk14400" + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# PHYSICS # +#.......................................# +incflo.gravity = 0. 0. -9.81 # Gravitational force (3D) +incflo.density = 1.0 # Reference density + +incflo.use_godunov = 1 +incflo.godunov_type = weno_z +incflo.diffusion_type = 1 +transport.viscosity = 1.0e-5 +transport.laminar_prandtl = 0.7 +transport.turbulent_prandtl = 0.3333 +turbulence.model = OneEqKsgsM84 + +incflo.physics = ABL Actuator +ICNS.source_terms = BoussinesqBuoyancy CoriolisForcing ABLMeanBoussinesq ActuatorForcing +TKE.source_terms = KsgsM84Src +BoussinesqBuoyancy.reference_temperature = 300.0 +CoriolisForcing.latitude = 41.3 +ABLForcing.abl_forcing_height = 90 +incflo.velocity = 6.928203230275509 4.0 0.0 + + +# Atmospheric boundary layer +ABL.temperature_heights = 0.0 700.0 800.0 1280.0 +ABL.temperature_values = 300.0 300.0 308.0 309.44 +ABL.reference_temperature = 300.0 +ABL.kappa = .40 +ABL.surface_roughness_z0 = 1.0E-4 +ABL.Uperiods = 25.0 +ABL.Vperiods = 25.0 +ABL.cutoff_height = 50.0 +ABL.deltaU = 1.0 +ABL.deltaV = 1.0 +ABL.normal_direction = 2 +ABL.perturb_ref_height = 50.0 +ABL.perturb_temperature = false +ABL.perturb_velocity = true +ABL.stats_output_format = netcdf +ABL.stats_output_frequency = 1 +ABL.surface_temp_flux = 0.00 +ABL.wall_shear_stress_type = "Moeng" + +ABL.bndry_file = "/projects/ssc/amr_precursors/b_abl_neutral_lowTI_redo/bndry_files" +ABL.bndry_io_mode = 1 +ABL.bndry_planes = ylo xlo # I'm (Paul) adding this but not sure if I have to +ABL.bndry_var_names = velocity temperature tke + + +# Output boundary files +ABL.bndry_planes = ylo xlo +ABL.bndry_output_start_time = 7200.0 +ABL.bndry_var_names = velocity temperature tke +ABL.bndry_output_format = native +ABL.stats_output_frequency = 1 +ABL.stats_output_format = netcdf + +# Whether to use helics +helics.activated = true +helics.broker_port = 32000 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# ADAPTIVE MESH REFINEMENT # +#.......................................# +amr.n_cell = 512 512 128 # Grid cells at coarsest AMRlevel +amr.max_level = 0 # Max AMR level in hierarchy + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# GEOMETRY # +#.......................................# +geometry.prob_lo = 0. 0. 0. # Lo corner coordinates +geometry.prob_hi = 5120. 5120. 1280. # Hi corner coordinates +geometry.is_periodic = 0 0 0 +xlo.type = "mass_inflow" +xlo.density = 1.0 +xlo.temperature = 0.0 # value required but ignored +xlo.tke = 0.0 +xhi.type = "pressure_outflow" + +ylo.type = "mass_inflow" +ylo.density = 1.0 +ylo.temperature = 0.0 +ylo.tke = 0.0 +yhi.type = "pressure_outflow" + +# Boundary conditions +zlo.type = "wall_model" +zlo.tke_type = "zero_gradient" + +zhi.type = "slip_wall" +zhi.temperature_type = "fixed_gradient" +zhi.temperature = 0.003 # tracer is used to specify potential temperature gradient + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# VERBOSITY # +#.......................................# +incflo.verbose = 0 # incflo_level + + + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# SAMPLING # +#.......................................# +incflo.post_processing = samplingPlane samplingLine + +samplingPlane.output_frequency = 600 +samplingPlane.labels = z_plane +samplingPlane.fields = velocity temperature +samplingPlane.z_plane.type = PlaneSampler +samplingPlane.z_plane.axis1 = 5110 0.0 0.0 +samplingPlane.z_plane.axis2 = 0.0 5110 0.0 +samplingPlane.z_plane.origin = 5.0 5.0 0.0 +samplingPlane.z_plane.num_points = 512 512 +samplingPlane.z_plane.normal = 0.0 0.0 1.0 +samplingPlane.z_plane.offsets = 5.0 85.0 155.0 255.0 + + +samplingLine.output_frequency = 1 +samplingLine.labels = z_line +samplingLine.fields = velocity temperature +samplingLine.z_line.type = LineSampler +samplingLine.z_line.num_points = 128 +samplingLine.z_line.start = 5.0 5.0 5.0 +samplingLine.z_line.end = 5.0 5.0 1275.0 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# TURBINES # +#.......................................# + +Actuator.type = UniformCtDisk +Actuator.UniformCtDisk.rotor_diameter = 126.0 +Actuator.UniformCtDisk.hub_height = 90.0 +Actuator.UniformCtDisk.thrust_coeff = 0.0 0.0 1.132034888 0.999470963 0.917697381 0.860849503 0.815371198 0.811614904 0.807939328 0.80443352 0.800993851 0.79768116 0.794529244 0.791495834 0.788560434 0.787217182 0.787127977 0.785839257 0.783812219 0.783568108 0.783328285 0.781194418 0.777292539 0.773464375 0.769690236 0.766001924 0.762348072 0.758760824 0.755242872 0.751792927 0.748434131 0.745113997 0.717806682 0.672204789 0.63831272 0.610176496 0.585456847 0.563222111 0.542912273 0.399312061 0.310517829 0.248633226 0.203543725 0.169616419 0.143478955 0.122938861 0.106515296 0.093026095 0.081648606 0.072197368 0.064388275 0.057782745 0.0 0.0 +Actuator.UniformCtDisk.wind_speed = 0.0 2.9 3.0 4.0 5.0 6.0 7.0 7.1 7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9 8.0 9.0 10.0 10.1 10.2 10.3 10.4 10.5 10.6 10.7 10.8 10.9 11.0 11.1 11.2 11.3 11.4 11.5 11.6 11.7 11.8 11.9 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 25.1 50.0 +Actuator.UniformCtDisk.epsilon = 10.0 +Actuator.UniformCtDisk.density = 1.225 +Actuator.UniformCtDisk.diameters_to_sample = 1.0 +Actuator.UniformCtDisk.num_points_r = 20 +Actuator.UniformCtDisk.num_points_t = 5 + + +Actuator.labels = T00 T01 +# T02 T03 T04 T05 T06 T07 +Actuator.JoukowskyDisk.thrust_coeff = 8.1672e-01 7.9044e-01 7.8393e-01 7.8624e-01 7.8824e-01 7.8942e-01 7.8902e-01 7.8740e-01 7.8503e-01 7.8237e-01 7.7955e-01 7.7583e-01 7.7583e-01 7.7583e-01 7.7583e-01 7.7583e-01 7.7583e-01 7.7583e-01 7.7583e-01 7.7583e-01 7.7583e-01 7.7583e-01 7.7583e-01 7.7583e-01 7.7583e-01 7.7583e-01 7.7583e-01 7.7583e-01 7.6922e-01 7.4270e-01 5.5949e-01 4.6163e-01 3.8786e-01 3.2901e-01 2.8093e-01 2.4114e-01 2.0795e-01 1.8010e-01 1.5663e-01 1.3679e-01 1.1995e-01 1.0562e-01 9.3384e-02 8.2908e-02 7.3910e-02 6.6159e-02 5.9463e-02 5.3662e-02 4.8622e-02 4.4230e-02 +Actuator.JoukowskyDisk.wind_speed = 3.0000e+00 3.5495e+00 4.0679e+00 4.5539e+00 5.0064e+00 5.4244e+00 5.8069e+00 6.1530e+00 6.4619e+00 6.7330e+00 6.9655e+00 7.1589e+00 7.3128e+00 7.4269e+00 7.5009e+00 7.5345e+00 7.5412e+00 7.5883e+00 7.6757e+00 7.8031e+00 7.9702e+00 8.1767e+00 8.4221e+00 8.7059e+00 9.0273e+00 9.3856e+00 9.7800e+00 1.0210e+01 1.0659e+01 1.0673e+01 1.1170e+01 1.1699e+01 1.2259e+01 1.2848e+01 1.3465e+01 1.4109e+01 1.4778e+01 1.5471e+01 1.6185e+01 1.6921e+01 1.7674e+01 1.8445e+01 1.9231e+01 2.0030e+01 2.0841e+01 2.1661e+01 2.2489e+01 2.3323e+01 2.4160e+01 2.5000e+01 +Actuator.JoukowskyDisk.rpm = 5.0000e+00 5.0000e+00 5.0000e+00 5.0000e+00 5.0000e+00 5.0000e+00 5.0000e+00 5.0000e+00 5.0000e+00 5.0000e+00 5.0000e+00 5.0861e+00 5.1954e+00 5.2765e+00 5.3290e+00 5.3529e+00 5.3577e+00 5.3912e+00 5.4532e+00 5.5437e+00 5.6625e+00 5.8092e+00 5.9836e+00 6.1851e+00 6.4135e+00 6.6681e+00 6.9483e+00 7.2535e+00 7.4992e+00 7.4992e+00 7.4992e+00 7.4992e+00 7.4992e+00 7.4992e+00 7.4992e+00 7.4992e+00 7.4992e+00 7.4992e+00 7.4992e+00 7.4992e+00 7.4992e+00 7.4992e+00 7.4992e+00 7.4992e+00 7.4992e+00 7.4992e+00 7.4992e+00 7.4992e+00 7.4992e+00 7.4992e+00 +Actuator.JoukowskyDisk.rotor_diameter = 240.0 +Actuator.JoukowskyDisk.hub_height = 150.0 +Actuator.JoukowskyDisk.output_frequency = 100 +Actuator.JoukowskyDisk.diameters_to_sample = 1.0 +Actuator.JoukowskyDisk.num_points_r = 40 +Actuator.JoukowskyDisk.num_points_t = 5 +Actuator.JoukowskyDisk.num_blades = 3 +Actuator.JoukowskyDisk.use_tip_correction = true +Actuator.JoukowskyDisk.use_root_correction = true +Actuator.JoukowskyDisk.epsilon = 5.0 +Actuator.JoukowskyDisk.vortex_core_size = 24.0 + +Actuator.UniformCtDisk.yaw = 240.0 + +Actuator.T00.base_position = 0.0 0.0 0.0 +Actuator.T01.base_position = 1000.0 0.0 0.0 diff --git a/tests/test_inputs/amr_standin_data.csv b/tests/test_inputs/amr_standin_data.csv new file mode 100644 index 00000000..306ff2f8 --- /dev/null +++ b/tests/test_inputs/amr_standin_data.csv @@ -0,0 +1,11 @@ +,time,amr_wind_speed,amr_wind_direction, +0,0.0,8.0,240.0 +1,1.0,10.0,240.0 +2,2.0,10.0,240.0 +3,3.0,8.0,240.0 +4,4.0,8.0,270.0 +5,5.0,8.0,270.0 +6,6.0,8.0,270.0 +7,7.0,8.0,270.0 +8,8.0,10.0,270.0 +9,9.0,10.0,270.0 \ No newline at end of file