diff --git a/openfast_toolbox/fastfarm/AMRWindSimulation.py b/openfast_toolbox/fastfarm/AMRWindSimulation.py index 67e3edf..d899f2a 100644 --- a/openfast_toolbox/fastfarm/AMRWindSimulation.py +++ b/openfast_toolbox/fastfarm/AMRWindSimulation.py @@ -23,16 +23,43 @@ def __init__(self, wts:dict, buffer_hr = 0.6, ds_hr = None, ds_lr = None, dt_hr = None, dt_lr = None, - mod_wake = None): + mod_wake = None, + level_lr = 0, + level_hr = -1): ''' Values from the AMR-Wind input file Inputs: - * dt: this should be a fixed dt value from the LES run - * incflo_velocity_hh: velocity vector, specifically at hub height - * buffer_lr: buffer for [xmin, xmax, ymin, ymax, zmax] in low-res box, in D - * buffer_hr: buffer for all directions (constant) in high-res box, in D - * mod_wake: Wake formulations within FAST.Farm. 1:polar; 2:curl; 3:cartesian. + ------ + dt: scalar + This should be a fixed dt value from the LES run + prob_lo, prob_hi: tuple + Extents of the AMR-Wind domain + n_cell: tuple of integers + Number of cells from AMR-Wind input file + max_level: int + Max level on the AMR-Wind grid. Used for high-res box placement + incflo_velocity_hh: tuple + Velocity vector, specifically at hub height + buffer_lr: 5-position vector + Buffer for [xmin, xmax, ymin, ymax, zmax] in low-res box, in D + buffer_hr: scalar + Buffer for all directions (constant) in high-res box, in D + ds_hr, ds_lr: scalar + Spatial resolution of the high-res and low-res boxes. No need to set + ahead of time, but optional to fix when you know the resolution you need + dt_hr, dt_lr: scalar + Same as above, but temporal resolution + mod_wake: int + Wake formulations within FAST.Farm. 1:polar; 2:curl; 3:cartesian. + level_lr: int + Grid level used for low-resolution box placement. Useful if staggered + grid approach is used for background mesh + level_hr: int + Grid level used for high-resolution boxes placement. Defaults to + highest level available + ''' + # Process inputs self.wts = wts self.dt = dt @@ -50,6 +77,8 @@ def __init__(self, wts:dict, self.dt_hr = dt_hr self.dt_lr = dt_lr self.mod_wake = mod_wake + self.level_lr = level_lr + self.level_hr = level_hr # Placeholder variables, to be calculated by FFCaseCreation self.output_frequency_lr = None @@ -73,9 +102,13 @@ def __init__(self, wts:dict, self._calc_simple_params() self._calc_sampling_params() + # Get execution time + from datetime import datetime + self.curr_datetime = datetime.now().strftime("%Y-%m-%d %H:%M:%S") + def __repr__(self): - s = f'<{type(self).__name__} object>\n' + s = f'<{type(self).__name__} object>\n\n' s += f'Requested parameters:\n' s += f' - Wake model: {self.mod_wake} (1:Polar; 2:Curl; 3:Cartesian)\n' s += f' - Extent of high-res boxes: {self.extent_high} D to each side\n' @@ -84,12 +117,13 @@ def __repr__(self): s += f'\n' s += f'LES parameters:\n' s += f' - velocity hub height: {self.incflo_velocity_hh} m/s\n' - s += f' - ds LES: ({self.dx0}, {self.dy0}, {self.dz0}) m\n' + s += f' - ds LES at low-res level (level {self.level_lr}): ({self.dx_at_lr_level}, {self.dy_at_lr_level}, {self.dz_at_lr_level}) m\n' + s += f' - ds LES at high-res level (level {self.level_hr}): ({self.dx_at_hr_level}, {self.dy_at_hr_level}, {self.dz_at_hr_level}) m\n' s += f' - dt LES: {self.dt} s\n' s += f' - Extents: ({self.prob_hi[0]-self.prob_lo[0]}, {self.prob_hi[1]-self.prob_lo[1]}, {self.prob_hi[2]-self.prob_lo[2]}) m\n' - s += f' - x: {self.prob_lo[0]}:{self.dx0}:{self.prob_hi[0]} m,\t ({self.n_cell[0]} points)\n' - s += f' - y: {self.prob_lo[1]}:{self.dy0}:{self.prob_hi[1]} m,\t ({self.n_cell[1]} points)\n' - s += f' - z: {self.prob_lo[2]}:{self.dz0}:{self.prob_hi[2]} m,\t ({self.n_cell[2]} points)\n' + s += f' - x: {self.prob_lo[0]}:{self.dx0}:{self.prob_hi[0]} m,\t ({self.n_cell[0]} points at level 0)\n' + s += f' - y: {self.prob_lo[1]}:{self.dy0}:{self.prob_hi[1]} m,\t ({self.n_cell[1]} points at level 0)\n' + s += f' - z: {self.prob_lo[2]}:{self.dz0}:{self.prob_hi[2]} m,\t ({self.n_cell[2]} points at level 0)\n' s += f'\n' s += f'Low-res domain: \n' @@ -97,9 +131,9 @@ def __repr__(self): s += f' - dt low: {self.dt_low_les} s (with LES dt = {self.dt} s, output frequency is {self.output_frequency_lr})\n' s += f' - Sampling labels: {self.sampling_labels_lr}\n' s += f' - Extents: ({self.xdist_lr}, {self.ydist_lr}, {self.zdist_lr}) m\n' - s += f' - x: {self.xlow_lr}:{self.ds_low_les}:{self.xhigh_lr} m,\t ({self.nx_lr} points)\n' - s += f' - y: {self.ylow_lr}:{self.ds_low_les}:{self.yhigh_lr} m,\t ({self.ny_lr} points)\n' - s += f' - z: {self.zlow_lr}:{self.ds_low_les}:{self.zhigh_lr} m,\t ({self.nz_lr} points)\n' + s += f' - x: {self.xlow_lr}:{self.ds_low_les}:{self.xhigh_lr} m,\t ({self.nx_lr} points at level {self.level_lr})\n' + s += f' - y: {self.ylow_lr}:{self.ds_low_les}:{self.yhigh_lr} m,\t ({self.ny_lr} points at level {self.level_lr})\n' + s += f' - z: {self.zlow_lr}:{self.ds_low_les}:{self.zhigh_lr} m,\t ({self.nz_lr} points at level {self.level_lr})\n' s += f'\n' s += f'High-res domain: \n' @@ -107,11 +141,11 @@ def __repr__(self): s += f' - dt high: {self.dt_high_les} s (with LES dt = {self.dt} s, output frequency is {self.output_frequency_hr})\n' s += f' - Sampling labels: {self.sampling_labels_hr}\n' for t in np.arange(len(self.hr_domains)): - s += f" - Turbine {t}\n" + s += f" - Turbine {t}, base located at ({self.wts[t]['x']:.2f}, {self.wts[t]['y']:.2f}, {self.wts[t]['z']:.2f}), hub height of {self.wts[t]['zhub']} m\n" s += f" - Extents: ({self.hr_domains[t]['xdist_hr']}, {self.hr_domains[t]['ydist_hr']}, {self.hr_domains[t]['zdist_hr']}) m\n" - s += f" - x: {self.hr_domains[t]['xlow_hr']}:{self.ds_high_les}:{self.hr_domains[t]['xhigh_hr']} m,\t ({self.hr_domains[t]['nx_hr']} points)\n" - s += f" - y: {self.hr_domains[t]['ylow_hr']}:{self.ds_high_les}:{self.hr_domains[t]['yhigh_hr']} m,\t ({self.hr_domains[t]['ny_hr']} points)\n" - s += f" - z: {self.hr_domains[t]['zlow_hr']}:{self.ds_high_les}:{self.hr_domains[t]['zhigh_hr']} m,\t ({self.hr_domains[t]['nz_hr']} points)\n" + s += f" - x: {self.hr_domains[t]['xlow_hr']}:{self.ds_high_les}:{self.hr_domains[t]['xhigh_hr']} m,\t ({self.hr_domains[t]['nx_hr']} points at level {self.level_hr})\n" + s += f" - y: {self.hr_domains[t]['ylow_hr']}:{self.ds_high_les}:{self.hr_domains[t]['yhigh_hr']} m,\t ({self.hr_domains[t]['ny_hr']} points at level {self.level_hr})\n" + s += f" - z: {self.hr_domains[t]['zlow_hr']}:{self.ds_high_les}:{self.hr_domains[t]['zhigh_hr']} m,\t ({self.hr_domains[t]['nz_hr']} points at level {self.level_hr})\n" return s @@ -134,6 +168,28 @@ def _checkInputs(self): if self.mod_wake not in [1,2,3]: raise ValueError (f'mod_wake parameter can only be 1 (polar), 2 (curl), or 3 (cartesian). Received {self.mod_wake}.') + if self.level_hr == -1: + self.level_hr = self.max_level + + # fmax and cmax from input farm should be >0 + # Check that level_lr is >= 0 + # check that level_hr is <=self.max_level + + + # Flags of given/calculated spatial resolution for warning/error printing purposes + self.given_ds_hr = False + self.given_ds_lr = False + warn_msg = "" + if self.ds_hr is not None: + warn_msg += f"--- WARNING: HIGH-RES SPATIAL RESOLUTION GIVEN. CONVERTING FATAL ERRORS ON HIGH-RES BOXES CHECKS TO WARNINGS. ---\n" + self.given_ds_hr = True + if self.ds_lr is not None: + warn_msg += f"--- WARNING: LOW-RES SPATIAL RESOLUTION GIVEN. CONVERTING FATAL ERRORS ON LOW-RES BOX CHECKS TO WARNINGS. ---\n" + self.given_ds_lr = True + print(f'{warn_msg}\n') + a=1 + + def _calc_simple_params(self): ''' Calculate simulation parameters, given only AMR-Wind inputs @@ -142,17 +198,30 @@ def _calc_simple_params(self): self.dx0 = (self.prob_hi[0] - self.prob_lo[0]) / self.n_cell[0] self.dy0 = (self.prob_hi[1] - self.prob_lo[1]) / self.n_cell[1] self.dz0 = (self.prob_hi[2] - self.prob_lo[2]) / self.n_cell[2] - self.ds0_max = max(self.dx0, self.dy0, self.dz0) - # Grid resolution at finest refinement level - self.dx_refine = self.dx0/(2**self.max_level) - self.dy_refine = self.dy0/(2**self.max_level) - self.dz_refine = self.dz0/(2**self.max_level) - self.ds_refine_max = max(self.dx_refine, self.dy_refine, self.dz_refine) + # Create grid resolutions at every level. Each array position corresponds to one level + # Access using res[1]['dx'], where `1` is the desired level + res = {} + for lvl in np.arange(0,self.max_level+1): + res[lvl] = {'dx':self.dx0/(2**lvl), 'dy':self.dy0/(2**lvl), 'dz':self.dz0/(2**lvl)} + + # Get grid resolutions at the levels of each box + self.dx_at_lr_level = res[self.level_lr]['dx'] + self.dy_at_lr_level = res[self.level_lr]['dy'] + self.dz_at_lr_level = res[self.level_lr]['dz'] + + self.dx_at_hr_level = res[self.level_hr]['dx'] + self.dy_at_hr_level = res[self.level_hr]['dy'] + self.dz_at_hr_level = res[self.level_hr]['dz'] + + # Get the maximum ds at each boxes' level + self.ds_max_at_lr_level = max(res[self.level_lr]['dx'], res[self.level_lr]['dy'], res[self.level_lr]['dz']) + self.ds_max_at_hr_level = max(res[self.level_hr]['dx'], res[self.level_hr]['dy'], res[self.level_hr]['dz']) # Hub height wind speed self.vhub = np.sqrt(self.incflo_velocity_hh[0]**2 + self.incflo_velocity_hh[1]**2) + def _calc_sampling_params(self): ''' Calculate parameters for sampling planes @@ -163,6 +232,7 @@ def _calc_sampling_params(self): self._calc_grid_resolution() self._calc_grid_placement() + def _calc_sampling_labels(self): ''' Calculate labels for AMR-Wind sampling @@ -180,42 +250,38 @@ def _calc_sampling_labels(self): self.sampling_labels_hr = sampling_labels_hr + def _calc_sampling_time(self): ''' Calculate timestep values and AMR-Wind plane sampling frequency ''' - ## High resolution domain, dt_high_les - fmax_max = 0 - for turbkey in self.wts: - fmax_max = max(0, self.wts[turbkey]['fmax']) + # Get some paramters from farm input + self.fmax_max = max(turb['fmax'] for turb in self.wts.values()) + self.cmeander_min = min(turb['Cmeander'] for turb in self.wts.values()) + self.Dwake_min = min(turb['D'] for turb in self.wts.values()) # Approximate D_wake as D_rotor + self.cmax_min = min(turb['cmax'] for turb in self.wts.values()) + + + # High resolution domain, dt_high_les if self.dt_hr is None: # Calculate dt of high-res per guidelines - dt_hr_max = 1 / (2 * fmax_max) + dt_hr_max = 1 / (2 * self.fmax_max) self.dt_high_les = getMultipleOf(dt_hr_max, multipleof=self.dt) # Ensure dt_hr is a multiple of the AMR-Wind timestep else: - # dt of high-res is given + # The dt of high-res is given self.dt_high_les = self.dt_hr if self.dt_high_les < self.dt: raise ValueError(f"AMR-Wind timestep {self.dt} too coarse for high resolution domain! AMR-Wind timestep must be at least {self.dt_high_les} sec.") - ## Low resolution domain, dt_low_les - cmeander_min = float("inf") - Dwake_min = float("inf") - for turbkey in self.wts: - cmeander_min = min(cmeander_min, self.wts[turbkey]['Cmeander']) - Dwake_min = min(Dwake_min, self.wts[turbkey]['D']) # Approximate D_wake as D_rotor - cmax_min = float("inf") - for turbkey in self.wts: - self.cmax_min = min(cmax_min, self.wts[turbkey]['cmax']) - - if self.mod_wake == 1: + # Low resolution domain, dt_low_les + if self.mod_wake == 1: # Polar self.dr = self.cmax_min - dt_lr_max = cmeander_min * Dwake_min / (10 * self.vhub) + dt_lr_max = self.cmeander_min * self.Dwake_min / (10 * self.vhub) else: # mod_wake == 2 or 3 - self.dr = Dwake_min/10 + self.dr = self.Dwake_min/15 dt_lr_max = self.dr / (2* self.vhub) @@ -232,72 +298,109 @@ def _calc_sampling_time(self): if self.dt_high_les > self.dt_low_les: raise ValueError(f"Low resolution timestep ({self.dt_low_les}) is finer than high resolution timestep ({self.dt_high_les})!") - ## Sampling frequency + + # Sampling frequency self.output_frequency_hr = int(np.floor(round(self.dt_high_les/self.dt,4))) self.output_frequency_lr = getMultipleOf(self.dt_low_les/self.dt, multipleof=self.output_frequency_hr) if self.output_frequency_lr % self.output_frequency_hr != 0: raise ValueError(f"Low resolution output frequency of {self.output_frequency_lr} not a multiple of the high resolution frequency {self.output_frequency_hr}!") + + def _calc_grid_resolution(self): ''' Calculate sampling grid resolutions + + Uses level_lr and level_hr to determine which AMR-Wind level is to be used to sample these boxes ''' - ## High resolution domain, ds_hr - # ASSUME: FAST.Farm HR zone lies within the region of maxmum AMR-Wind grid refinement - # NOTE: ds_hr is calculated independent of any x/y/z requirements, - # just blade chord length requirements + # Calculate ds_hr and ds_lr if not specified as inputs + if self.ds_hr is None: + self.ds_hr = self._calc_grid_resolution_hr() + + if self.ds_lr is None: + self.ds_lr = self._calc_grid_resolution_lr() + + # Set the computed or speficied spatial resolutions + self.ds_high_les = self.ds_hr + self.ds_low_les = self.ds_lr + + + # Perform some checks + if self.ds_high_les < self.ds_max_at_hr_level: + raise ValueError(f"AMR-Wind grid spacing of {self.ds_max_at_hr_level} m at the high-res box level of {self.level_hr} is too coarse for "\ + f"the high resolution domain! AMR-Wind grid spacing at level {self.level_hr} must be at least {self.ds_high_les} m.") + + if self.ds_low_les < self.ds_max_at_lr_level: + error_msg = f"AMR-Wind grid spacing of {self.ds_max_at_lr_level} at the low-res box level of {self.level_lr} is too coarse for "\ + f"the low resolution domain! AMR-Wind grid spacing at level {self.level_lr} must be at least {self.ds_low_les} m. "\ + f"If you can afford to have {self.ds_low_les} m on AMR-Wind for the low-res box, do so. If you cannot, add `ds_lr={self.ds_max_at_lr_level}` " \ + f"to the call to `AMRWindSimulation`. Note that sampled values will no longer be at the cell centers, as you will be requesting "\ + f"sampling at {self.ds_low_les} m while the underlying grid will be at {self.ds_max_at_lr_level} m.\n --- SUPRESSING FURTHER ERRORS ---" + if self.given_ds_hr: + print(f'WARNING: {error_msg}') + else: + raise ValueError(error_msg) + + if self.ds_low_les % self.ds_high_les != 0: + raise ValueError(f"Low resolution grid spacing of {self.ds_low_les} m not a multiple of the high resolution grid spacing {self.ds_high_les} m.") + + + + def _calc_grid_resolution_hr(self): + ds_hr_max = self.cmax_min - if self.ds_hr is None: # Calculate ds_hr if it is not specified as an input - if ds_hr_max < self.ds_refine_max: - raise ValueError(f"AMR-Wind grid spacing of {self.ds_refine_max} is too coarse for high resolution domain! The high-resolution domain requires "\ - f"AMR-Wind grid spacing to be at least {ds_hr_max} m. If a coarser high-res domain is acceptable, then manually specify the "\ - f"high-resolution grid spacing to be at least {self.ds_refine_max} with ds_hr = {self.ds_refine_max}.") - self.ds_high_les = getMultipleOf(ds_hr_max, multipleof=self.ds_refine_max) # Ensure that ds_hr is a multiple of the refined AMR-Wind grid spacing + if ds_hr_max < self.ds_max_at_hr_level: + raise ValueError(f"AMR-Wind grid spacing of {self.ds_max_at_hr_level} m at the high-res box level of {self.level_hr} is too coarse for "\ + f"high resolution domain! The high-resolution domain requires AMR-Wind grid spacing to be at least {ds_hr_max} m. If a "\ + f"coarser high-res domain is acceptable, then manually specify the high-resolution grid spacing to be at least "\ + f"{self.ds_max_at_hr_level} with ds_hr = {self.ds_max_at_hr_level}.") + ds_high_les = getMultipleOf(ds_hr_max, multipleof=self.ds_max_at_hr_level) # Ensure that ds_hr is a multiple of the refined AMR-Wind grid spacing - self.ds_hr = self.ds_high_les - else: - self.ds_high_les = self.ds_hr + return ds_high_les + + + def _calc_grid_resolution_lr(self): + + ds_lr_max = self.cmeander_min * self.Dwake_min * self.vhub / 150 + # The expression above is the same as + # For polar wake model: ds_lr_max = self.dt_low_les * self.vhub**2 / 15 + # For curled wake model: ds_lr_max = self.cmeander_max * self.dt_low_les * self.vhub**2 / 5 + + ds_low_les = getMultipleOf(ds_lr_max, multipleof=self.ds_hr) + print(f"Low-res spatial resolution should be at least {ds_lr_max:.2f} m, but since it needs to be a multiple of high-res "\ + f"resolution of {self.ds_hr}, we pick ds_low to be {ds_low_les} m") + + #self.ds_lr = self.ds_low_les + return ds_low_les - if self.ds_high_les < self.ds_refine_max: - raise ValueError(f"AMR-Wind fine grid spacing {self.ds_refine_max} too coarse for high resolution domain! AMR-Wind grid spacing must be at least {self.ds_high_les} m.") - ## Low resolution domain, ds_lr (s = x/y/z) - # ASSUME: FAST.Farm LR zone uses Level 0 AMR-Wind grid spacing - # NOTE: ds_lr is calculated independent of any x/y/z requirements, - # just time step and velocity requiements - if self.ds_lr is None: - ds_lr_max = self.dt_low_les * self.vhub**2 / 15 - self.ds_low_les = getMultipleOf(ds_lr_max, multipleof=self.ds_hr) # ds_hr is already a multiple of the AMR-Wind grid spacing, so here we need to make sure ds_lr is a multiple of ds_hr - self.ds_lr = self.ds_low_les - else: - self.ds_low_les = self.ds_lr - if self.ds_low_les < self.ds0_max: - raise ValueError(f"AMR-Wind coarse grid spacing {self.ds0_max} too coarse for high resolution domain! AMR-Wind grid spacing must be at least {self.ds_low_les} m.") - if self.ds_low_les % self.ds_high_les != 0: - raise ValueError(f"Low resolution grid spacing of {self.ds_low_les} not a multiple of the high resolution grid spacing {self.ds_high_les}!") def _calc_grid_placement(self): ''' Calculate placement of sampling grids ''' - self._calc_hr_grid_placement() - self._calc_lr_grid_placement() + self._calc_grid_placement_hr() + self._calc_grid_placement_lr() - def _calc_hr_grid_placement(self): + + def _calc_grid_placement_hr(self): ''' - Calculate placement of high resolution grids + Calculate placement of high resolution grids based on desired level level_hr + ''' - ### ~~~~~~~~~ Calculate high resolution grid placement ~~~~~~~~~ + + # Calculate high resolution grid placement hr_domains = {} for turbkey in self.wts: wt_x = self.wts[turbkey]['x'] wt_y = self.wts[turbkey]['y'] - wt_z = self.wts[turbkey]['zhub'] + wt_z = self.wts[turbkey]['z'] + wt_h = self.wts[turbkey]['zhub'] wt_D = self.wts[turbkey]['D'] # Calculate minimum/maximum HR domain extents @@ -305,7 +408,8 @@ def _calc_hr_grid_placement(self): xhigh_hr_max = wt_x + self.buffer_hr * wt_D ylow_hr_min = wt_y - self.buffer_hr * wt_D yhigh_hr_max = wt_y + self.buffer_hr * wt_D - zhigh_hr_max = wt_z + self.buffer_hr * wt_D + zlow_hr_min = wt_z + zhigh_hr_max = wt_z + wt_h + self.buffer_hr * wt_D # Calculate the minimum/maximum HR domain coordinate lengths & number of grid cells xdist_hr_min = xhigh_hr_max - xlow_hr_min # Minumum possible length of x-extent of HR domain @@ -316,32 +420,44 @@ def _calc_hr_grid_placement(self): ydist_hr = self.ds_high_les * np.ceil(ydist_hr_min/self.ds_high_les) ny_hr = int(ydist_hr/self.ds_high_les) + 1 - zdist_hr = self.ds_high_les * np.ceil(zhigh_hr_max/self.ds_high_les) + zdist_hr_min = zhigh_hr_max - zlow_hr_min + zdist_hr = self.ds_high_les * np.ceil(zdist_hr_min/self.ds_high_les) nz_hr = int(zdist_hr/self.ds_high_les) + 1 # Calculate actual HR domain extent # NOTE: Sampling planes should measure at AMR-Wind cell centers, not cell edges - xlow_hr = getMultipleOf(xlow_hr_min, multipleof=self.ds_high_les) - 0.5*self.dx_refine + self.prob_lo[0]%self.ds_high_les + xlow_hr = getMultipleOf(xlow_hr_min, multipleof=self.ds_high_les) - 0.5*self.dx_at_hr_level+ self.prob_lo[0]%self.ds_high_les xhigh_hr = xlow_hr + xdist_hr - ylow_hr = getMultipleOf(ylow_hr_min, multipleof=self.ds_high_les) - 0.5*self.dy_refine + self.prob_lo[1]%self.ds_high_les + + ylow_hr = getMultipleOf(ylow_hr_min, multipleof=self.ds_high_les) - 0.5*self.dy_at_hr_level + self.prob_lo[1]%self.ds_high_les yhigh_hr = ylow_hr + ydist_hr - zlow_hr = 0.5 * self.dz0 / (2**self.max_level) + + #zlow_hr = zlow_hr_min + 0.5 * self.dz_at_hr_level + # !!!! 2024-07-23: changed to positive the 0.5*gridres below so that z>0 in flat terrain. The wfip files will be off + zlow_hr = getMultipleOf(zlow_hr_min, multipleof=self.ds_high_les) + 0.5*self.dz_at_hr_level + self.prob_lo[2]%self.ds_high_les zhigh_hr = zlow_hr + zdist_hr zoffsets_hr = np.arange(zlow_hr, zhigh_hr+self.ds_high_les, self.ds_high_les) - zlow_hr + # Check domain extents if xhigh_hr > self.prob_hi[0]: - raise ValueError(f"HR domain point {xhigh_hr} extends beyond maximum AMR-Wind x-extent!") + raise ValueError(f"Turbine {turbkey}: HR domain point {xhigh_hr} extends beyond maximum AMR-Wind x-extent!") + #print(f"Turbine {turbkey}: ERROR: HR domain point {xhigh_hr} extends beyond maximum AMR-Wind x-extent!") if xlow_hr < self.prob_lo[0]: - raise ValueError(f"HR domain point {xlow_hr} extends beyond minimum AMR-Wind x-extent!") + raise ValueError(f"Turbine {turbkey}: HR domain point {xlow_hr} extends beyond minimum AMR-Wind x-extent!") + #print(f"Turbine {turbkey}: ERROR: HR domain point {xlow_hr} extends beyond minimum AMR-Wind x-extent!") if yhigh_hr > self.prob_hi[1]: - raise ValueError(f"HR domain point {yhigh_hr} extends beyond maximum AMR-Wind y-extent!") + raise ValueError(f"Turbine {turbkey}: HR domain point {yhigh_hr} extends beyond maximum AMR-Wind y-extent!") + #print(f"Turbine {turbkey}: ERROR: HR domain point {yhigh_hr} extends beyond maximum AMR-Wind y-extent!") if ylow_hr < self.prob_lo[1]: - raise ValueError(f"HR domain point {ylow_hr} extends beyond minimum AMR-Wind y-extent!") + raise ValueError(f"Turbine {turbkey}: HR domain point {ylow_hr} extends beyond minimum AMR-Wind y-extent!") + #print(f"Turbine {turbkey}: ERROR: HR domain point {ylow_hr} extends beyond minimum AMR-Wind y-extent!") if zhigh_hr > self.prob_hi[2]: - raise ValueError(f"HR domain point {zhigh_hr} extends beyond maximum AMR-Wind z-extent!") + raise ValueError(f"Turbine {turbkey}: HR domain point {zhigh_hr} extends beyond maximum AMR-Wind z-extent!") + #print(f"Turbine {turbkey}: ERROR: HR domain point {zhigh_hr} extends beyond maximum AMR-Wind z-extent!") if zlow_hr < self.prob_lo[2]: - raise ValueError(f"HR domain point {zlow_hr} extends beyond minimum AMR-Wind z-extent!") + raise ValueError(f"Turbine {turbkey}: HR domain point {zlow_hr} extends beyond minimum AMR-Wind z-extent!") + #print(f"Turbine {turbkey}: ERROR: HR domain point {zlow_hr} extends beyond minimum AMR-Wind z-extent!") # Save out info for FFCaseCreation self.extent_high = self.buffer_hr*2 @@ -354,45 +470,44 @@ def _calc_hr_grid_placement(self): hr_domains[turbkey] = hr_turb_info self.hr_domains = hr_domains - def _calc_lr_grid_placement(self): + + def _calc_grid_placement_lr(self): ''' - Calculate placement of low resolution grid + Calculate placement of low resolution grid based on desired level level_lr ''' ### ~~~~~~~~~ Calculate low resolution grid placement ~~~~~~~~~ # Calculate minimum/maximum LR domain extents - wt_all_x_min = float("inf") # Minimum x-value of any turbine - wt_all_x_max = -1*float("inf") - wt_all_y_min = float("inf") - wt_all_y_max = -1*float("inf") - wt_all_z_max = -1*float("inf") # Tallest rotor disk point of any turbine - Drot_max = -1*float("inf") - for turbkey in self.wts: - wt_all_x_min = min(wt_all_x_min, self.wts[turbkey]['x']) - wt_all_x_max = max(wt_all_x_max, self.wts[turbkey]['x']) - wt_all_y_min = min(wt_all_y_min, self.wts[turbkey]['y']) - wt_all_y_max = max(wt_all_y_min, self.wts[turbkey]['y']) - wt_all_z_max = max(wt_all_z_max, self.wts[turbkey]['zhub'] + 0.5*self.wts[turbkey]['D']) - Drot_max = max(Drot_max, self.wts[turbkey]['D']) + wt_all_x_min = min(turb['x'] for turb in self.wts.values()) + wt_all_x_max = max(turb['x'] for turb in self.wts.values()) + wt_all_y_min = min(turb['y'] for turb in self.wts.values()) + wt_all_y_max = max(turb['y'] for turb in self.wts.values()) + wt_all_z_min = min(turb['z'] for turb in self.wts.values()) + wt_all_z_max = max(turb['z']+turb['zhub']+0.5*turb['D'] for turb in self.wts.values()) + D_max = max(turb['D'] for turb in self.wts.values()) - xlow_lr_min = wt_all_x_min - self.buffer_lr[0] * Drot_max - xhigh_lr_max = wt_all_x_max + self.buffer_lr[1] * Drot_max - ylow_lr_min = wt_all_y_min - self.buffer_lr[2] * Drot_max - yhigh_lr_max = wt_all_y_max + self.buffer_lr[3] * Drot_max - zhigh_lr_max = wt_all_z_max + self.buffer_lr[4] * Drot_max + xlow_lr_min = wt_all_x_min - self.buffer_lr[0] * D_max + xhigh_lr_max = wt_all_x_max + self.buffer_lr[1] * D_max + ylow_lr_min = wt_all_y_min - self.buffer_lr[2] * D_max + yhigh_lr_max = wt_all_y_max + self.buffer_lr[3] * D_max + zlow_lr_min = wt_all_z_min - 0.5 * D_max + zhigh_lr_max = wt_all_z_max + self.buffer_lr[4] * D_max + + # Carve out exception for flat terrain + if wt_all_z_min == 0: + zlow_lr_min = wt_all_z_min # Calculate the minimum/maximum LR domain coordinate lengths & number of grid cells xdist_lr_min = xhigh_lr_max - xlow_lr_min # Minumum possible length of x-extent of LR domain - self.xdist_lr = self.ds_low_les * np.ceil(xdist_lr_min/self.ds_low_les) # The `+ ds_lr` comes from the +1 to NS_LOW in Sec. 4.2.15.6.4.1.1 - # TODO: adjust xdist_lr calculation by also using `inflow_deg` - self.nx_lr = int(self.xdist_lr/self.ds_low_les) + 1 - ydist_lr_min = yhigh_lr_max - ylow_lr_min + zdist_lr_min = zhigh_lr_max - zlow_lr_min + + self.xdist_lr = self.ds_low_les * np.ceil(xdist_lr_min/self.ds_low_les) # The `+ ds_lr` comes from the +1 to NS_LOW in Sec. 4.2.15.6.4.1.1 self.ydist_lr = self.ds_low_les * np.ceil(ydist_lr_min/self.ds_low_les) - # TODO: adjust ydist_lr calculation by also using `inflow_deg` - self.ny_lr = int(self.ydist_lr/self.ds_low_les) + 1 + self.zdist_lr = self.ds_low_les * np.ceil(zdist_lr_min/self.ds_low_les) - self.zdist_lr = self.ds_low_les * np.ceil(zhigh_lr_max/self.ds_low_les) + self.nx_lr = int(self.xdist_lr/self.ds_low_les) + 1 # TODO: adjust xdist_lr calculation by also using `inflow_deg` + self.ny_lr = int(self.ydist_lr/self.ds_low_les) + 1 # TODO: adjust ydist_lr calculation by also using `inflow_deg` self.nz_lr = int(self.zdist_lr/self.ds_low_les) + 1 ## Calculate actual LR domain extent @@ -400,15 +515,16 @@ def _calc_lr_grid_placement(self): # NOTE: Should we use dx/dy/dz values here or ds_lr? # - AR: I think it's correct to use ds_lr to get to the xlow values, # but then offset by 0.5*amr_dx0 if need be - self.xlow_lr = getMultipleOf(xlow_lr_min, multipleof=self.ds_low_les) - 0.5*self.dx0 + self.prob_lo[0]%self.ds_low_les + self.xlow_lr = getMultipleOf(xlow_lr_min, multipleof=self.ds_low_les) - 0.5*self.dx_at_lr_level + self.prob_lo[0]%self.ds_low_les + self.ylow_lr = getMultipleOf(ylow_lr_min, multipleof=self.ds_low_les) - 0.5*self.dy_at_lr_level + self.prob_lo[1]%self.ds_low_les + self.zlow_lr = getMultipleOf(zlow_lr_min, multipleof=self.ds_low_les) + 0.5*self.dz_at_lr_level + self.prob_lo[2]%self.ds_low_les + self.xhigh_lr = self.xlow_lr + self.xdist_lr - self.ylow_lr = getMultipleOf(ylow_lr_min, multipleof=self.ds_low_les) - 0.5*self.dy0 + self.prob_lo[1]%self.ds_low_les self.yhigh_lr = self.ylow_lr + self.ydist_lr - self.zlow_lr = 0.5 * self.dz0 # Lowest z point is half the height of the lowest grid cell self.zhigh_lr = self.zlow_lr + self.zdist_lr self.zoffsets_lr = np.arange(self.zlow_lr, self.zhigh_lr+self.ds_low_les, self.ds_low_les) - self.zlow_lr - ## Check domain extents + # Check domain extents if self.xhigh_lr > self.prob_hi[0]: raise ValueError(f"LR domain point {self.xhigh_lr} extends beyond maximum AMR-Wind x-extent!") if self.xlow_lr < self.prob_lo[0]: @@ -422,54 +538,50 @@ def _calc_lr_grid_placement(self): if self.zlow_lr < self.prob_lo[2]: raise ValueError(f"LR domain point {self.zlow_lr} extends beyond minimum AMR-Wind z-extent!") - ## Check grid placement + # Check grid placement self._check_grid_placement() - ## Save out info for FFCaseCreation + # Save out info for FFCaseCreation self.extent_low = self.buffer_lr + def _check_grid_placement(self): ''' Check the values of parameters that were calculated by _calc_sampling_params ''' - ## Check that sampling grids are at cell centers - # Low resolution grid - amr_xgrid_level0 = np.arange(self.prob_lo[0], self.prob_hi[0], self.dx0) - amr_ygrid_level0 = np.arange(self.prob_lo[1], self.prob_hi[1], self.dy0) - amr_zgrid_level0 = np.arange(self.prob_lo[2], self.prob_hi[2], self.dz0) + # Calculate parameters of the low-res grid to allow checking + amr_xgrid_at_lr_level = np.arange(self.prob_lo[0], self.prob_hi[0], self.dx_at_lr_level) + amr_ygrid_at_lr_level = np.arange(self.prob_lo[1], self.prob_hi[1], self.dy_at_lr_level) + amr_zgrid_at_lr_level = np.arange(self.prob_lo[2], self.prob_hi[2], self.dz_at_lr_level) + + amr_xgrid_at_lr_level_cc = amr_xgrid_at_lr_level + 0.5*self.dx_at_lr_level # Cell-centered AMR-Wind x-grid + amr_ygrid_at_lr_level_cc = amr_ygrid_at_lr_level + 0.5*self.dy_at_lr_level + amr_zgrid_at_lr_level_cc = amr_zgrid_at_lr_level + 0.5*self.dz_at_lr_level - amr_xgrid_level0_cc = amr_xgrid_level0 + 0.5*self.dx0 # Cell-centered AMR-Wind x-grid - amr_ygrid_level0_cc = amr_ygrid_level0 + 0.5*self.dy0 - amr_zgrid_level0_cc = amr_zgrid_level0 + 0.5*self.dz0 - sampling_xgrid_lr = self.xlow_lr + self.ds_lr*np.arange(self.nx_lr) sampling_ygrid_lr = self.ylow_lr + self.ds_lr*np.arange(self.ny_lr) sampling_zgrid_lr = self.zlow_lr + self.zoffsets_lr - # TODO: These for loops could be replaced with a faster operation - for coord in sampling_xgrid_lr: - if coord not in amr_xgrid_level0_cc: - raise ValueError("Low resolution x-sampling grid is not cell cenetered with AMR-Wind's grid!") - for coord in sampling_ygrid_lr: - if coord not in amr_ygrid_level0_cc: - raise ValueError("Low resolution y-sampling grid is not cell cenetered with AMR-Wind's grid!") - for coord in sampling_zgrid_lr: - if coord not in amr_zgrid_level0_cc: - raise ValueError("Low resolution z-sampling grid is not cell cenetered with AMR-Wind's grid!") + # Check the low-res grid placement + self._check_grid_placement_single(sampling_xgrid_lr, amr_xgrid_at_lr_level_cc, boxstr='Low', dirstr='x') + self._check_grid_placement_single(sampling_ygrid_lr, amr_ygrid_at_lr_level_cc, boxstr='Low', dirstr='y') + self._check_grid_placement_single(sampling_zgrid_lr, amr_zgrid_at_lr_level_cc, boxstr='Low', dirstr='z') + + # High resolution grids (span the entire domain to make this check easier) - amr_xgrid_refine = np.arange(self.prob_lo[0], self.prob_hi[0], self.dx_refine) - amr_ygrid_refine = np.arange(self.prob_lo[1], self.prob_hi[1], self.dy_refine) - amr_zgrid_refine = np.arange(self.prob_lo[2], self.prob_hi[2], self.dz_refine) + amr_xgrid_at_hr_level = np.arange(self.prob_lo[0], self.prob_hi[0], self.dx_at_hr_level) + amr_ygrid_at_hr_level = np.arange(self.prob_lo[1], self.prob_hi[1], self.dy_at_hr_level) + amr_zgrid_at_hr_level = np.arange(self.prob_lo[2], self.prob_hi[2], self.dz_at_hr_level) - amr_xgrid_refine_cc = amr_xgrid_refine + 0.5*self.dx_refine - amr_ygrid_refine_cc = amr_ygrid_refine + 0.5*self.dy_refine - amr_zgrid_refine_cc = amr_zgrid_refine + 0.5*self.dz_refine + amr_xgrid_at_hr_level_cc = amr_xgrid_at_hr_level + 0.5*self.dx_at_hr_level + amr_ygrid_at_hr_level_cc = amr_ygrid_at_hr_level + 0.5*self.dy_at_hr_level + amr_zgrid_at_hr_level_cc = amr_zgrid_at_hr_level + 0.5*self.dz_at_hr_level for turbkey in self.hr_domains: - nx_hr = self.hr_domains[turbkey]['nx_hr'] - ny_hr = self.hr_domains[turbkey]['ny_hr'] + nx_hr = self.hr_domains[turbkey]['nx_hr'] + ny_hr = self.hr_domains[turbkey]['ny_hr'] xlow_hr = self.hr_domains[turbkey]['xlow_hr'] ylow_hr = self.hr_domains[turbkey]['ylow_hr'] @@ -477,32 +589,72 @@ def _check_grid_placement(self): sampling_ygrid_hr = ylow_hr + self.ds_hr*np.arange(ny_hr) sampling_zgrid_hr = self.hr_domains[turbkey]['zlow_hr'] + self.hr_domains[turbkey]['zoffsets_hr'] - # TODO: These for loops could be replaced with a faster operation - for coord in sampling_xgrid_hr: - if coord not in amr_xgrid_refine_cc: - raise ValueError("High resolution x-sampling grid is not cell cenetered with AMR-Wind's grid!") - for coord in sampling_ygrid_hr: - if coord not in amr_ygrid_refine_cc: - raise ValueError("High resolution y-sampling grid is not cell cenetered with AMR-Wind's grid!") - for coord in sampling_zgrid_hr: - if coord not in amr_zgrid_refine_cc: - raise ValueError("High resolution z-sampling grid is not cell cenetered with AMR-Wind's grid!") + # Check the high-res grid placement + self._check_grid_placement_single(sampling_xgrid_hr, amr_xgrid_at_hr_level_cc, boxstr='High', dirstr='x') + self._check_grid_placement_single(sampling_ygrid_hr, amr_ygrid_at_hr_level_cc, boxstr='High', dirstr='y') + self._check_grid_placement_single(sampling_zgrid_hr, amr_zgrid_at_hr_level_cc, boxstr='High', dirstr='z') + + def _check_grid_placement_single(self, sampling_xyzgrid_lhr, amr_xyzgrid_at_lhr_level_cc, boxstr, dirstr): + ''' + Generic function to check placement of x, y, z grids in low, high-res boxes + + Inputs + ------ + sampling_xyzgrid_lhr: np array + Actual sampling locations of {x,y,z}grid on either low or high res (`lhr`) + amr_xyzgrid_at_lhr_level_cc: np array + Actual AMR-Wind grid cell-center location of {x,y,z}grid at either low- or high-res levels + boxstr: str + Either 'Low', or 'High' + distr: str + Either 'x', 'y', or 'z' - def write_sampling_params(self, outdir=None): + ''' + + # Check if all values in sampling grid sampling_{x,y,z}grid_lr are part of AMR-Wind cell-centered values amr_{x,y,z}grid_at_lr_level_cc + is_sampling_xyzgrid_subset = set(sampling_xyzgrid_lhr).issubset(set(amr_xyzgrid_at_lhr_level_cc)) + + if is_sampling_xyzgrid_subset is False: + amr_index = np.argmin(np.abs(amr_xyzgrid_at_lhr_level_cc - sampling_xyzgrid_lhr[0])) + error_msg = f"{boxstr} resolution {dirstr}-sampling grid is not cell-centered with AMR-Wind's grid. \n "\ + f"{dirstr}-sampling grid: {sampling_xyzgrid_lhr[0]}, {sampling_xyzgrid_lhr[1]}, "\ + f"{sampling_xyzgrid_lhr[2]}, {sampling_xyzgrid_lhr[3]}, ... \n "\ + f"AMR-Wind grid (subset): {amr_xyzgrid_at_lhr_level_cc[amr_index ]}, {amr_xyzgrid_at_lhr_level_cc[amr_index+1]}, "\ + f"{amr_xyzgrid_at_lhr_level_cc[amr_index+2]}, {amr_xyzgrid_at_lhr_level_cc[amr_index+3]}, ..." + if self.given_ds_lr: + print(f'WARNING: {error_msg}') + else: + raise ValueError(error_msg) + + + + def write_sampling_params(self, out=None, overwrite=False): ''' Write out text that can be used for the sampling planes in an AMR-Wind input file - outdir: str - Input file to be written. If None, result is written to screen + out: str + Output path or full filename for Input file to be written + to. If None, result is written to screen. + overwrite: bool + If saving to a file, whether or not to overwrite potentially + existing file ''' + # Write time step information for consistenty with sampling frequency + s = f"time.fixed_dt = {self.dt}\n\n" + # Write flow velocity info for consistency + s += f"incflo.velocity = {self.incflo_velocity_hh[0]} {self.incflo_velocity_hh[1]} {self.incflo_velocity_hh[2]}\n\n\n" + # Write high-level info for sampling sampling_labels_lr_str = " ".join(str(item) for item in self.sampling_labels_lr) sampling_labels_hr_str = " ".join(str(item) for item in self.sampling_labels_hr) - s = f"# Sampling info generated by AMRWindSamplingCreation.py\n" + s += f"#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#\n" + s += f"# POST-Processing #\n" + s += f"#.......................................#\n" + s += f"# Sampling info generated by AMRWindSamplingCreation.py on {self.curr_datetime}\n" s += f"incflo.post_processing = {self.postproc_name_lr} {self.postproc_name_hr} # averaging\n\n\n" s += f"# ---- Low-res sampling parameters ----\n" @@ -533,6 +685,8 @@ def write_sampling_params(self, outdir=None): for turbkey in self.hr_domains: wt_x = self.wts[turbkey]['x'] wt_y = self.wts[turbkey]['y'] + wt_z = self.wts[turbkey]['z'] + wt_h = self.wts[turbkey]['zhub'] wt_D = self.wts[turbkey]['D'] if 'name' in self.wts[turbkey].keys(): wt_name = self.wts[turbkey]['name'] @@ -549,7 +703,7 @@ def write_sampling_params(self, outdir=None): zoffsets_hr = self.hr_domains[turbkey]['zoffsets_hr'] zoffsets_hr_str = " ".join(str(item) for item in zoffsets_hr) - s += f"\n# Turbine {wt_name} at (x,y) = ({wt_x}, {wt_y}), with D = {wt_D}, grid spacing = {self.ds_hr} m\n" + s += f"\n# Turbine {wt_name} with base at (x,y,z) = ({wt_x:.4f}, {wt_y:.4f}, {wt_z:.4f}), with hh = {wt_h}, D = {wt_D}, grid spacing = {self.ds_hr} m\n" s += f"{self.postproc_name_hr}.{sampling_name}.type = PlaneSampler\n" s += f"{self.postproc_name_hr}.{sampling_name}.num_points = {nx_hr} {ny_hr}\n" s += f"{self.postproc_name_hr}.{sampling_name}.origin = {xlow_hr:.4f} {ylow_hr:.4f} {zlow_hr:.4f}\n" # Round the float output @@ -559,18 +713,20 @@ def write_sampling_params(self, outdir=None): s += f"{self.postproc_name_hr}.{sampling_name}.offsets = {zoffsets_hr_str}\n" - if outdir is None: + if out is None: print(s) - else: - outfile = os.path.join(outdir, 'sampling_config.i') - if not os.path.exists(outdir): - print(f'Path {outdir} does not exist. Creating it') - os.makedirs(outdir) - if os.path.isfile(outfile): - raise FileExistsError(f"{str(outfile)} already exists! Aborting...") - - with open(outfile,"w") as out: - out.write(s) + return + elif os.path.isdir(out): + outfile = os.path.join(out, 'sampling_config.i') + else: + # full file given + outfile = out + if not overwrite: + if os.path.isfile(outfile): + raise FileExistsError(f"{str(outfile)} already exists! Aborting...") + + with open(outfile,"w") as out: + out.write(s) diff --git a/openfast_toolbox/fastfarm/FASTFarmCaseCreation.py b/openfast_toolbox/fastfarm/FASTFarmCaseCreation.py index 79da86e..7205583 100644 --- a/openfast_toolbox/fastfarm/FASTFarmCaseCreation.py +++ b/openfast_toolbox/fastfarm/FASTFarmCaseCreation.py @@ -12,6 +12,8 @@ def cosd(t): return np.cos(np.deg2rad(t)) def sind(t): return np.sin(np.deg2rad(t)) def checkIfExists(f): + if os.path.basename(f) == 'unused': + return True if os.path.isfile(f): return True else: @@ -82,7 +84,6 @@ def __init__(self, nSeeds = 6, seedValues = None, LESpath = None, - sweepWakeSteering = False, sweepYawMisalignment = False, refTurb_rot = 0, verbose = 0): @@ -142,8 +143,6 @@ def __init__(self, Full path of the LES data, if driven by LES. If None, the setup will be for TurbSim inflow. LESpath can be a single path, or a list of paths of the same length as the sweep in conditions. For example, if TIvalue=[8,10,12], then LESpath can be 3 paths, related to each condition. - sweepWakeSteering: bool - Whether or not to perform a sweep with wake steering sweepYawMisalignment: bool Whether or not to perform a sweep with and without yaw misalignment perturbations refTurb_rot: int @@ -175,12 +174,10 @@ def __init__(self, self.EDmodel = EDmodel self.nSeeds = nSeeds self.LESpath = LESpath - self.sweepWS = sweepWakeSteering self.sweepYM = sweepYawMisalignment self.seedValues = seedValues self.refTurb_rot = refTurb_rot self.verbose = verbose - self.mod_wake = mod_wake self.attempt = 1 @@ -335,6 +332,8 @@ def _checkInputs(self): if self.yaw_init is None: yaw = np.ones((1,self.nTurbines))*0 self.yaw_init = np.repeat(yaw, len(self.inflow_deg), axis=0) + if not isinstance(self.yaw_init, np.ndarray): + self.yaw_init = np.array(self.yaw_init) # Check TI values if given in percent for t in self.TIvalue: @@ -376,10 +375,24 @@ def _checkInputs(self): if len(self.inflow_deg) != len(self.yaw_init): raise ValueError(f'One row for each inflow angle should be given in yaw_init. '\ f'Currently {len(self.inflow_deg)} inflow angle(s) and {len(self.yaw_init)} yaw entrie(s)') - if self.ADmodel is None: + + # Check reduced-order models + if self.ADmodel is None or self.ADmodel == 'ADyn': self.ADmodel = np.tile(['ADyn'],(1,self.nTurbines)) - if self.EDmodel is None: + elif self.ADmodel == 'ADsk': + self.ADmodel = np.tile(['ADsk'],(1,self.nTurbines)) + else: + raise ValueError(f"Aerodynamic model {self.ADmodel} not recognized. Options are None or 'ADyn' "\ + f"for AeroDyn and 'ADsk' for AeroDisk") + + if self.EDmodel is None or self.EDmodel == 'FED' or self.EDmodel == 'ED': self.EDmodel = np.tile(['FED'],(1,self.nTurbines)) + elif self.EDmodel == 'SED': + self.EDmodel = np.tile(['SED'],(1,self.nTurbines)) + else: + raise ValueError(f"Elastic model {self.EDmodel} not recognized. Options are None, 'FED', or 'ED' "\ + f"for ElastoDyn and 'SED' for Simplified ElastoDyn") + if np.shape(self.ADmodel) != np.shape(self.EDmodel): raise ValueError('Every case should have the aerodynamic and elastic model selected. The number of cases '\ '(lines) in `ADmodel` and `EDmodel` should be the same') @@ -405,9 +418,11 @@ def _checkInputs(self): # Check LES parameters if self.LESpath is None: self.inflowStr = 'TurbSim' + self.Mod_AmbWind = 3 else: if isinstance(self.LESpath,str): self.LESpath = [self.LESpath]*len(self.vhub) self.inflowStr = 'LES' + self.Mod_AmbWind = 1 for p in self.LESpath: if not os.path.isdir(p): raise ValueError (f'The LES path {p} does not exist') @@ -439,7 +454,6 @@ def _checkInputs(self): raise ValueError(f'The grid resolution dS_High should not be greater than dS_Low on the LES side') - # Check the reference turbine for rotation if self.refTurb_rot >= self.nTurbines: raise ValueError(f'The index for the reference turbine for the farm to be rotated around is greater than the number of turbines') @@ -448,6 +462,7 @@ def _checkInputs(self): self.templateFilesCreatedBool = False self.TSlowBoxFilesCreatedBool = False self.TShighBoxFilesCreatedBool = False + self.hasController = False @@ -510,7 +525,6 @@ def _create_dir_structure(self): for case in range(self.nCases): # Recover information about current case for directory naming purposes inflow_deg_ = self.allCases['inflow_deg' ].sel(case=case).values - #wakeSteering_ = self.allCases['wakeSteering' ].sel(case=case).values misalignment_ = self.allCases['misalignment' ].sel(case=case).values nADyn_ = self.allCases['nFullAeroDyn' ].sel(case=case).values nFED_ = self.allCases['nFulllElastoDyn'].sel(case=case).values @@ -520,8 +534,6 @@ def _create_dir_structure(self): ndigits = len(str(self.nCases)) caseStr = f"Case{case:0{ndigits}d}_wdir{f'{int(inflow_deg_):+03d}'.replace('+','p').replace('-','m')}" # Add standard sweeps to the case name - if self.sweepWS: - caseStr += f"_WS{str(wakeSteering_).lower()}" if self.sweepYM: caseStr += f"_YM{str(misalignment_).lower()}" if self.sweepEDmodel: @@ -584,7 +596,7 @@ def copyTurbineFilesForEachCase(self, writeFiles=True): self.HydroDynFile.write(os.path.join(currPath, self.HDfilename)) # Write updated DISCON - if writeFiles: + if writeFiles and self.hasController: shutilcopy2_untilSuccessful(os.path.join(self.templatePath,self.controllerInputfilename), os.path.join(currPath,self.controllerInputfilename)) @@ -593,28 +605,42 @@ def copyTurbineFilesForEachCase(self, writeFiles=True): # the compilation process and it is likely that a very long string with the full path will get cut. So we need to # give the relative path. We give the path as the current one, so here we create a link to ensure it will work # regardless of how the controller was compiled. There is no harm in having this extra link even if it's not needed. - notepath = os.getcwd(); os.chdir(self.path) - for seed in range(self.nSeeds): - try: - src = os.path.join('../', self.controllerInputfilename) - dst = os.path.join(self.condDirList[cond], self.caseDirList[case], f'Seed_{seed}', self.controllerInputfilename) - if writeFiles: - os.symlink(src, dst) - except FileExistsError: - pass - os.chdir(notepath) + if self.hasController: + notepath = os.getcwd(); os.chdir(self.path) + for seed in range(self.nSeeds): + try: + src = os.path.join('..', self.controllerInputfilename) + dst = os.path.join(self.condDirList[cond], self.caseDirList[case], f'Seed_{seed}', self.controllerInputfilename) + if writeFiles: + os.symlink(src, dst) + except FileExistsError: + pass + os.chdir(notepath) - # Write InflowWind files. For FAST.FARM, the IW file needs to be inside the Seed* directories. If running standalone openfast, - # it needs to be on the same level as the fst file. Here, we copy to both places so that the workflow is general + # Write InflowWind files. For FAST.FARM with TS, the IW file needs to be inside the Seed* directories. If running standalone openfast + # or Mod_AmbWind==1 (VTK), it needs to be on the same level as the fst file. Here, we copy to both places so that the workflow is general self.InflowWindFile['WindType'] = 3 self.InflowWindFile['PropagationDir'] = 0 self.InflowWindFile['Filename_BTS'] = '"./TurbSim"' + self.InflowWindFile['NWindVel'] = 1 + self.InflowWindFile['WindVxiList'] = 0 # Sampling relative to the local reference frame + self.InflowWindFile['WindVyiList'] = 0 + self.InflowWindFile['WindVziList'] = self.allCases.sel(case=case, turbine=0)['zhub'].values + if writeFiles: self.InflowWindFile.write( os.path.join(currPath,self.IWfilename)) - for seed in range(self.nSeeds): - self.InflowWindFile.write( os.path.join(currPath,f'Seed_{seed}',self.IWfilename)) + if self.Mod_AmbWind == 3: # only for TS-driven cases + for seed in range(self.nSeeds): + self.InflowWindFile.write( os.path.join(currPath,f'Seed_{seed}',self.IWfilename)) + # Before starting the loop, print once the info about the controller is no controller is present + if not self.hasController: + if self.verbose>=1: + print(f" No controller given through libdiscon/DLL. ", + f"Using `VSContrl` {self.ServoDynFile['VSContrl']} from the template files.") + + # Loop through all turbines of current condition and case for t in range(self.nTurbines): # Recover info about the current turbine in CondXX_*/CaseYY_ yaw_deg_ = self.allCases.sel(case=case, turbine=t)['yaw'].values @@ -657,23 +683,44 @@ def copyTurbineFilesForEachCase(self, writeFiles=True): # Update each turbine's ServoDyn self.ServoDynFile['YawNeut'] = yaw_deg_ + yaw_mis_deg_ - self.ServoDynFile['DLL_FileName'] = f'"{self.DLLfilepath}{t+1}.so"' + if self.hasController: + self.ServoDynFile['VSContrl'] = 5 + self.ServoDynFile['DLL_FileName'] = f'"{self.DLLfilepath}{t+1}.so"' + self.ServoDynFile['DLL_InFile'] = f'"{self.controllerInputfilename}"' + else: + self.ServoDynFile['DLL_FileName'] = f'"unused"' + self.ServoDynFile['DLL_InFile'] = f'"unused"' if writeFiles: self.ServoDynFile.write( os.path.join(currPath,f'{self.SrvDfilename}{t+1}_mod.dat')) # Update each turbine's OpenFAST input self.turbineFile['TMax'] = self.tmax self.turbineFile['CompInflow'] = 1 # 1: InflowWind; 2: OpenFoam (fully coupled; not VTK input to FF) - + + if self.SubDfilename == 'unused': + self.turbineFile['CompSub'] = 0 + else: + self.turbineFile['CompSub'] = 1 + + if self.HDfilename == 'unused': + self.turbineFile['CompHydro'] = 0 + else: + self.turbineFile['CompHydro'] = 1 + + #if self.MDfilename == 'unused': + # self.turbineFile['CompMooring'] = 0 + #else: + # self.turbineFile['CompMooring'] = 1 + if EDmodel_ == 'FED': self.turbineFile['CompElast'] = 1 # 1: full ElastoDyn; 2: full ElastoDyn + BeamDyn; 3: Simplified ElastoDyn - self.turbineFile['CompSub'] = 1 - self.turbineFile['CompHydro'] = 1 + #self.turbineFile['CompSub'] = 1 + #self.turbineFile['CompHydro'] = 1 self.turbineFile['EDFile'] = f'"./{self.EDfilename}{t+1}_mod.dat"' elif EDmodel_ == 'SED': self.turbineFile['CompElast'] = 3 # 1: full ElastoDyn; 2: full ElastoDyn + BeamDyn; 3: Simplified ElastoDyn - self.turbineFile['CompSub'] = 0 # need to be disabled with SED - self.turbineFile['CompHydro'] = 0 # need to be disabled with SED + #self.turbineFile['CompSub'] = 0 # need to be disabled with SED + #self.turbineFile['CompHydro'] = 0 # need to be disabled with SED self.turbineFile['IntMethod'] = 3 self.turbineFile['EDFile'] = f'"./{self.SEDfilename}{t+1}_mod.dat"' self.turbineFile['BDBldFile(1)'] = f'"{self.BDfilepath}"' @@ -689,8 +736,8 @@ def copyTurbineFilesForEachCase(self, writeFiles=True): self.turbineFile['AeroFile'] = f'"{self.ADskfilepath}"' if writeFiles: if t==0: shutilcopy2_untilSuccessful(self.coeffTablefilepath, os.path.join(currPath,self.coeffTablefilename)) - self.turbineFile['ServoFile'] = f'"./{self.SrvDfilename}{t+1}_mod.dat"' - self.turbineFile['HydroFile'] = f'"./{self.HDfilename}"' + self.turbineFile['ServoFile'] = f'"{self.SrvDfilename}{t+1}_mod.dat"' + self.turbineFile['HydroFile'] = f'"{self.HDfilename}"' self.turbineFile['SubFile'] = f'"{self.SubDfilepath}"' self.turbineFile['MooringFile'] = f'"unused"' self.turbineFile['IceFile'] = f'"unused"' @@ -729,8 +776,11 @@ def _were_all_turbine_files_copied(self): _ = checkIfExists(os.path.join(currPath, self.HDfilename)) if not _: return False - _ = checkIfExists(os.path.join(currPath,self.controllerInputfilename)) - if not _: return False + + if self.hasController: + _ = checkIfExists(os.path.join(currPath,self.controllerInputfilename)) + if not _: return False + _ = checkIfExists( os.path.join(currPath,self.IWfilename)) if not _: return False @@ -809,7 +859,9 @@ def setTemplateFilename(self, self.BDfilepath = "unused"; self.BDfilename = "unused" self.bladefilename = "unused"; self.bladefilepath = "unused" self.towerfilename = "unused"; self.towerfilepath = "unused" - + self.libdisconfilepath = None + self.controllerInputfilename = None + self.coeffTablefilename = None if templatePath is None: print(f'--- WARNING: No template files given. Complete setup will not be possible') @@ -841,6 +893,12 @@ def checkIfExists(f): checkIfExists(self.SEDfilepath) self.SEDfilename = SEDfilename + if HDfilename is not None and HDfilename != 'unused': + raise ValueError (f'Simplified ElastoDyn is not compatible with HydroDyn. Set HDfilename to None. ') + if SubDfilename is not None and SubDfilename != 'unused': + raise ValueError (f'Simplified ElastoDyn is not compatible with SubDyn. Set SubDfilename to None. ') + + if HDfilename is not None and HDfilename != 'unused': if not HDfilename.endswith('.dat'): raise ValueError (f'The HydroDyn filename should end in `.dat`.') @@ -868,6 +926,14 @@ def checkIfExists(f): self.ADskfilepath = os.path.join(self.templatePath,ADskfilename) checkIfExists(self.ADskfilepath) self.ADskfilename = ADskfilename + self.hasController = False + + if coeffTablefilename is not None and coeffTablefilename != 'unused': + if not coeffTablefilename.endswith('.csv'): + raise ValueError (f'The performance table `coeffTablefilename` file is needed for AeroDisk and should end in "*.csv"') + self.coeffTablefilepath = os.path.join(templatePath, coeffTablefilename) + checkIfExists(self.coeffTablefilepath) + self.coeffTablefilename = coeffTablefilename if SubDfilename is not None and SubDfilename != 'unused': if not SubDfilename.endswith('.dat'): @@ -916,6 +982,7 @@ def checkIfExists(f): self.libdisconfilepath = libdisconfilepath checkIfExists(self.libdisconfilepath) self._create_copy_libdiscon() + self.hasController = True if controllerInputfilename is not None: if not controllerInputfilename.endswith('.IN'): @@ -923,13 +990,9 @@ def checkIfExists(f): self.controllerInputfilepath = os.path.join(self.templatePath, controllerInputfilename) checkIfExists(self.controllerInputfilepath) self.controllerInputfilename = controllerInputfilename - - if coeffTablefilename is not None and coeffTablefilename != 'unused': - if not coeffTablefilename.endswith('.csv'): - raise ValueError (f'The performance table `coeffTablefilename` file should end in "*.csv"') - self.coeffTablefilepath = os.path.join(templatePath, coeffTablefilename) - checkIfExists(self.coeffTablefilepath) - self.coeffTablefilename = coeffTablefilename + else: + if self.hasController: + raise ValueError (f'libdiscon has been given but no controller input filename. Stopping.') if turbsimLowfilepath is not None: if not turbsimLowfilepath.endswith('.inp'): @@ -1034,7 +1097,7 @@ def _create_all_cond(self): def _create_all_cases(self): - # Generate the different "cases" (inflow angle, and misalignment and wakesteer bools). + # Generate the different "cases" (inflow angle and yaw misalignment bools). # If misalignment true, then the actual yaw is yaw[turb]=np.random.uniform(low=-8.0, high=8.0). # Set sweep bools and multipliers @@ -1044,6 +1107,7 @@ def _create_all_cases(self): self.sweepEDmodel = False self.sweepADmodel = False else: + print(f"Sweeps in AD and ED enabled.") self.sweepEDmodel = True self.sweepADmodel = True @@ -1257,20 +1321,24 @@ def TS_low_slurm_prepare(self, slurmfilepath): memory_per_cpu = int(150000/self.nSeeds) # Change job name (for convenience only) - _ = subprocess.call(f"sed -i 's|#SBATCH --job-name=lowBox|#SBATCH --job-name=lowBox_{os.path.basename(self.path)}|g' {self.slurmfilename_low}", cwd=self.path, shell=True) + sed_command = f"sed -i 's|^#SBATCH --job-name=lowBox|#SBATCH --job-name=lowBox_{os.path.basename(self.path)}|g' {self.slurmfilename_low}" + _ = subprocess.call(sed_command, cwd=self.path, shell=True) # Change the path inside the script to the desired one - sed_command = f"sed -i 's|/projects/shellwind/rthedin/Task2_2regis|{self.path}|g' {self.slurmfilename_low}" + sed_command = f"""sed -i "s|^basepath.*|basepath='{self.path}'|g" {self.slurmfilename_low}""" _ = subprocess.call(sed_command, cwd=self.path, shell=True) # Change number of nodes values - _ = subprocess.call(f"sed -i 's|#SBATCH --nodes=2|#SBATCH --nodes={int(np.ceil(self.nConditions*self.nSeeds/6))}|g' {self.slurmfilename_low}", cwd=self.path, shell=True) + sed_command = f"sed -i 's|^#SBATCH --nodes.*|#SBATCH --nodes={int(np.ceil(self.nConditions*self.nSeeds/6))}|g' {self.slurmfilename_low}" + _ = subprocess.call(sed_command, cwd=self.path, shell=True) # Change memory per cpu - _ = subprocess.call(f"sed -i 's|--mem-per-cpu=25000M|--mem-per-cpu={memory_per_cpu}M|g' {self.slurmfilename_low}", cwd=self.path, shell=True) + sed_command = f"sed -i 's|--mem-per-cpu=25000M|--mem-per-cpu={memory_per_cpu}M|g' {self.slurmfilename_low}" + _ = subprocess.call(sed_command, cwd=self.path, shell=True) # Assemble list of conditions and write it listtoprint = "' '".join(self.condDirList) sed_command = f"""sed -i "s|^condList.*|condList=('{listtoprint}')|g" {self.slurmfilename_low}""" _ = subprocess.call(sed_command, cwd=self.path, shell=True) # Change the number of seeds - _ = subprocess.call(f"sed -i 's|nSeeds=6|nSeeds={self.nSeeds}|g' {self.slurmfilename_low}", cwd=self.path, shell=True) + sed_command = f"sed -i 's|^nSeeds.*|nSeeds={self.nSeeds}|g' {self.slurmfilename_low}" + _ = subprocess.call(sed_command, cwd=self.path, shell=True) if self.nSeeds > 6: @@ -1304,7 +1372,7 @@ def TS_low_createSymlinks(self): for case in range(self.nCases): for seed in range(self.nSeeds): try: - src = os.path.join('../../../..', self.condDirList[cond], f'Seed_{seed}', 'Low.bts') + src = os.path.join('..', '..', '..', '..', self.condDirList[cond], f'Seed_{seed}', 'Low.bts') dst = os.path.join(self.condDirList[cond], self.caseDirList[case], f'Seed_{seed}', 'TurbSim', 'Low.bts') os.symlink(src, dst) except FileExistsError: @@ -1319,7 +1387,7 @@ def getDomainParameters(self): if self.verbose>1: print(' Running a TurbSim setup once to get domain extents') self.TS_low_setup(writeFiles=False, runOnce=True) - # Figure out how many (and which) high boxes actually need to be executed. Remember that wake steering, yaw misalignment, SED/ADsk models, + # Figure out how many (and which) high boxes actually need to be executed. Remember that yaw misalignment, SED/ADsk models, # and sweep in yaw do not require extra TurbSim runs self.nHighBoxCases = len(np.unique(self.inflow_deg)) # some wind dir might be repeated for sweep on yaws @@ -1531,16 +1599,20 @@ def TS_high_slurm_prepare(self, slurmfilepath): shutil.copy2(slurmfilepath, os.path.join(self.path, self.slurmfilename_high)) # Change job name (for convenience only) - _ = subprocess.call(f"sed -i 's|#SBATCH --job-name=highBox|#SBATCH --job-name=highBox_{os.path.basename(self.path)}|g' {self.slurmfilename_high}", cwd=self.path, shell=True) + sed_command = f"sed -i 's|^#SBATCH --job-name.*|#SBATCH --job-name=highBox_{os.path.basename(self.path)}|g' {self.slurmfilename_high}" + _ = subprocess.call(sed_command, cwd=self.path, shell=True) # Change the path inside the script to the desired one - sed_command = f"sed -i 's|/projects/shellwind/rthedin/Task2_2regis|{self.path}|g' {self.slurmfilename_high}" + sed_command = f"""sed -i "s|^basepath.*|basepath='{self.path}'|g" {self.slurmfilename_high}""" _ = subprocess.call(sed_command, cwd=self.path, shell=True) # Change number of turbines - _ = subprocess.call(f"sed -i 's|nTurbines=12|nTurbines={self.nTurbines}|g' {self.slurmfilename_high}", cwd=self.path, shell=True) + sed_command = f"sed -i 's|^nTurbines.*|nTurbines={self.nTurbines}|g' {self.slurmfilename_high}" + _ = subprocess.call(sed_command, cwd=self.path, shell=True) # Change number of seeds - _ = subprocess.call(f"sed -i 's|nSeeds=6|nSeeds={self.nSeeds}|g' {self.slurmfilename_high}", cwd=self.path, shell=True) + set_command = f"sed -i 's|^nSeeds.*|nSeeds={self.nSeeds}|g' {self.slurmfilename_high}" + _ = subprocess.call(sed_command, cwd=self.path, shell=True) # Change number of nodes values - _ = subprocess.call(f"sed -i 's|#SBATCH --nodes=3|#SBATCH --nodes={int(np.ceil(ntasks/36))}|g' {self.slurmfilename_high}", cwd=self.path, shell=True) + sed_command = f"sed -i 's|^#SBATCH --nodes.*|#SBATCH --nodes={int(np.ceil(ntasks/36))}|g' {self.slurmfilename_high}" + _ = subprocess.call(sed_command, cwd=self.path, shell=True) # Assemble list of conditions and write it listtoprint = "' '".join(self.condDirList) sed_command = f"""sed -i "s|^condList.*|condList=('{listtoprint}')|g" {self.slurmfilename_high}""" @@ -1573,7 +1645,7 @@ def TS_high_slurm_submit(self, qos='normal', A=None, t=None, p=None): def TS_high_create_symlink(self): - # Create symlink of all the high boxes for the cases with wake steering and yaw misalignment. These are the "repeated" boxes + # Create symlink of all the high boxes for the cases with yaw misalignment. These are the "repeated" boxes if self.verbose>0: print(f'Creating symlinks for all the high-resolution boxes') @@ -1609,7 +1681,7 @@ def TS_high_create_symlink(self): # Now that we have the correct arrays, we perform the loop on the turbines and seeds for t in range(self.nTurbines): for seed in range(self.nSeeds): - src = os.path.join('../../../..', self.condDirList[cond], self.caseDirList[src_case], f'Seed_{seed}', 'TurbSim', f'HighT{t+1}.bts') + src = os.path.join('..', '..', '..', '..', self.condDirList[cond], self.caseDirList[src_case], f'Seed_{seed}', 'TurbSim', f'HighT{t+1}.bts') dst = os.path.join(self.condDirList[cond], self.caseDirList[case], f'Seed_{seed}', 'TurbSim', f'HighT{t+1}.bts') try: @@ -1780,8 +1852,8 @@ def _FF_setup_LES(self, seedsToKeep=1): ff_file = FASTInputFile(outputFSTF) # Open output file and change additional values manually or make sure we have the correct ones - ff_file['InflowFile'] = f'"./{self.IWfilename}"' - ff_file['Mod_AmbWind'] = 1 # 1: LES boxes; 2: single TurbSim; 3: multiple TurbSim + ff_file['InflowFile'] = f'"unused"' + ff_file['Mod_AmbWind'] = self.Mod_AmbWind # LES ff_file['TMax'] = self.tmax # LES-related parameters @@ -1880,7 +1952,7 @@ def _FF_setup_TS(self): ff_file = FASTInputFile(outputFSTF) ff_file['InflowFile'] = f'"./{self.IWfilename}"' #ff_file['DT']=1.0 - ff_file['Mod_AmbWind'] = 3 # 1: LES boxes; 2: single TurbSim; 3: multiple TurbSim + ff_file['Mod_AmbWind'] = self.Mod_AmbWind # 3: multiple TurbSim ff_file['TMax'] = self.tmax # Super controller @@ -2082,7 +2154,7 @@ def FF_slurm_prepare(self, slurmfilepath): sed_command = f"""sed -i "s|^fastfarmbin.*|fastfarmbin='{self.ffbin}'|g" {fname}""" _ = subprocess.call(sed_command, cwd=self.path, shell=True) # Change the path inside the script to the desired one - sed_command = f"sed -i 's|/projects/shellwind/rthedin/Task2_2regis|{self.path}|g' {fname}" + sed_command = f"""sed -i "s|^basepath.*|basepath='{self.path}'|g" {fname}""" _ = subprocess.call(sed_command, cwd=self.path, shell=True) # Write condition sed_command = f"""sed -i "s|^cond.*|cond='{self.condDirList[cond]}'|g" {fname}""" @@ -2156,7 +2228,7 @@ def FF_check_output(self): - def plot(self, figsize=(15,7), fontsize=14, saveFig=True, returnFig=False, figFormat='png'): + def plot(self, figsize=(14,7), fontsize=13, saveFig=True, returnFig=False, figFormat='png', showTurbNumber=False): import matplotlib.pyplot as plt fig, ax = plt.subplots(figsize=figsize) @@ -2181,7 +2253,7 @@ def plot(self, figsize=(15,7), fontsize=14, saveFig=True, returnFig=False, figFo for j, inflow in enumerate(self.wts_rot_ds['inflow_deg']): ax.set_prop_cycle(None) # Reset the colormap for every inflow for i, currTurbine in enumerate(self.wts_rot_ds.turbine): - color = next(ax._get_lines.prop_cycler)['color'] + color = ax._get_lines.get_next_color() dst = self.wts_rot_ds.sel(turbine=currTurbine, inflow_deg=inflow) @@ -2198,6 +2270,9 @@ def plot(self, figsize=(15,7), fontsize=14, saveFig=True, returnFig=False, figFo # plot turbine location ax.scatter(dst.x, dst.y, s=dst.D/6, c=color, marker='o') #, label=f'WT{i+1}') + if j==0 and showTurbNumber: + # Show turbine number + ax.annotate(f'T{str(currTurbine.values+1)}', (dst.x, dst.y), textcoords="offset points", xytext=(4,4), ha='center') # plot turbine disk accoding to all yaws in current wdir allyaw_currwdir = self.allCases.where(self.allCases['inflow_deg']==inflow,drop=True).sel(turbine=currTurbine)['yaw'] @@ -2210,6 +2285,7 @@ def plot(self, figsize=(15,7), fontsize=14, saveFig=True, returnFig=False, figFo # plot convex hull of farm (or line) for given inflow turbs = self.wts_rot_ds.sel(inflow_deg=inflow)[['x','y']].to_array().transpose() try: + from scipy.spatial import ConvexHull hull = ConvexHull(turbs) for simplex in hull.simplices: ax.plot(turbs[simplex, 0], turbs[simplex, 1], 'gray', alpha=alphas[j], label=f'Inflow {inflow.values} deg') @@ -2221,7 +2297,7 @@ def plot(self, figsize=(15,7), fontsize=14, saveFig=True, returnFig=False, figFo # Remove duplicate entries from legend handles, labels = plt.gca().get_legend_handles_labels() by_label = dict(zip(labels, handles)) - plt.legend(by_label.values(), by_label.keys(), loc='upper left', bbox_to_anchor=(1.02,1.015), fontsize=fontsize) + plt.legend(by_label.values(), by_label.keys(), loc='upper left', bbox_to_anchor=(1.02,1.015), fontsize=fontsize, ncols=int(self.nTurbines/25)) ax.set_xlabel("x [m]", fontsize=fontsize) ax.set_ylabel("y [m]", fontsize=fontsize) diff --git a/openfast_toolbox/fastfarm/examples/Ex1_TurbSimInputSetup.py b/openfast_toolbox/fastfarm/examples/Ex1_TurbSimInputSetup.py index 61d5f3e..f255e53 100644 --- a/openfast_toolbox/fastfarm/examples/Ex1_TurbSimInputSetup.py +++ b/openfast_toolbox/fastfarm/examples/Ex1_TurbSimInputSetup.py @@ -16,8 +16,8 @@ MyDir=os.path.dirname(__file__) # --- Define parameters necessary for this script -OldTSFile = os.path.join(MyDir, 'SampleFiles/TestCase.inp' ) # Template file used for TurbSim, need to exist -NewTSFile = os.path.join(MyDir, 'SampleFiles/_TestCase_mod.inp') # New file that will be written +OldTSFile = os.path.join(MyDir, 'SampleFiles', 'TestCase.inp' ) # Template file used for TurbSim, need to exist +NewTSFile = os.path.join(MyDir, 'SampleFiles', '_TestCase_mod.inp') # New file that will be written D = 77.0 # Turbine diameter (m) HubHt = 78.045 # Hub Height (m) Vhub = 6 # mean wind speed at hub height (m/s) diff --git a/openfast_toolbox/fastfarm/examples/Ex2_FFarmInputSetup.py b/openfast_toolbox/fastfarm/examples/Ex2_FFarmInputSetup.py index 96895c5..9e1a9f1 100644 --- a/openfast_toolbox/fastfarm/examples/Ex2_FFarmInputSetup.py +++ b/openfast_toolbox/fastfarm/examples/Ex2_FFarmInputSetup.py @@ -20,8 +20,8 @@ if __name__ == '__main__': # --- FAST Farm input files - templateFSTF = os.path.join(MyDir,'SampleFiles/TestCase.fstf') # template file used for FastFarm input file, need to exist - outputFSTF = os.path.join(MyDir,'SampleFiles/_TestCase_mod.fstf')# new file that will be written + templateFSTF = os.path.join(MyDir,'SampleFiles', 'TestCase.fstf') # template file used for FastFarm input file, need to exist + outputFSTF = os.path.join(MyDir,'SampleFiles', '_TestCase_mod.fstf')# new file that will be written # --- Parameters for TurbSim Extent D = 77.0 # Turbine diameter (m) hubHeight = 78.045 # Hub Height (m) @@ -29,7 +29,7 @@ extent_YZ_high = 1.2 # y-extent of high res box in diamter around turbine location chord_max = 5 # maximum blade chord (m). Turbine specific. Cmeander = 1.9 # Meandering constant (-) - BTSFilename = os.path.join(MyDir,'SampleFiles/TestCase.bts') # TurbSim Box to be used in FAST.Farm simulation, need to exist. + BTSFilename = os.path.join(MyDir,'SampleFiles', 'TestCase.bts') # TurbSim Box to be used in FAST.Farm simulation, need to exist. # --- Layout xWT = [0.0, 265.] # x positions of turbines yWT = [0.0, 50.0] # y postitions of turbines diff --git a/openfast_toolbox/fastfarm/examples/Ex3_FFarmCompleteSetup.py b/openfast_toolbox/fastfarm/examples/Ex3_FFarmCompleteSetup.py index f2903fc..f4914f6 100644 --- a/openfast_toolbox/fastfarm/examples/Ex3_FFarmCompleteSetup.py +++ b/openfast_toolbox/fastfarm/examples/Ex3_FFarmCompleteSetup.py @@ -60,7 +60,7 @@ def main(): # ----------- Turbine parameters # Set the yaw of each turbine for wind dir. One row for each wind direction. - yaw_init = [ [0,0,0,0,0,0,0,0,0,0,0,0] ] + yaw_init = None # ----------- Low- and high-res boxes parameters # Should match LES if comparisons are to be made; otherwise, set desired values @@ -87,7 +87,7 @@ def main(): # ----------- Template files templatePath = '/full/path/where/template/files/are' - # Put 'unused' to any input that is not applicable to your case + # Put None on any input that is not applicable to your case # Files should be in templatePath EDfilename = 'ElastoDyn.T' SEDfilename = 'SimplifiedElastoDyn.T' @@ -97,7 +97,7 @@ def main(): ADskfilename = 'AeroDisk.dat' SubDfilename = 'SubDyn.dat' IWfilename = 'InflowWind.dat' - BDfilepath = 'unused' + BDfilepath = None bladefilename = 'Blade.dat' towerfilename = 'Tower.dat' turbfilename = 'Model.T' @@ -125,8 +125,8 @@ def main(): case = FFCaseCreation(path, wts, tmax, zbot, vhub, shear, TIvalue, inflow_deg, dt_high_les, ds_high_les, extent_high, dt_low_les, ds_low_les, extent_low, - ffbin, mod_wake, yaw_init, - nSeeds=nSeeds, LESpath=LESpath, + ffbin=ffbin, mod_wake=mod_wake, yaw_init=yaw_init, + nSeeds=nSeeds, LESpath=LESpath, refTurb_rot=refTurb_rot, verbose=1) case.setTemplateFilename(templatePath, EDfilename, SEDfilename, HDfilename, SrvDfilename, ADfilename,