From d274f46e11f6e08f0656e9e5f1a3608b4185f2aa Mon Sep 17 00:00:00 2001 From: jefalon Date: Fri, 1 Mar 2024 16:07:51 -0700 Subject: [PATCH 01/15] added ability to integrate turbines locally --- .github/workflows/CI_every_PR.yml | 14 +- .github/workflows/CI_every_PR_parallel.yml | 50 +++++ tests/9-Regression/ALMx2_Unsteady.yaml | 2 +- windse/ProblemManager.py | 75 +++++-- windse/SolverManager.py | 12 +- windse/input_schema.yaml | 12 ++ windse/turbine_types/ActuatorDisk2D.py | 2 + windse/turbine_types/ActuatorLine.py | 4 +- windse/turbine_types/GenericTurbine.py | 2 +- windse/wind_farm_types/GenericWindFarm.py | 237 +++++++++++++++------ windse_driver/driver.py | 3 + 11 files changed, 313 insertions(+), 100 deletions(-) create mode 100644 .github/workflows/CI_every_PR_parallel.yml diff --git a/.github/workflows/CI_every_PR.yml b/.github/workflows/CI_every_PR.yml index 3bdb0768..15cbeb9a 100644 --- a/.github/workflows/CI_every_PR.yml +++ b/.github/workflows/CI_every_PR.yml @@ -41,13 +41,13 @@ jobs: # cd tests/9-Regression # mpirun -n 2 windse run ALM_Unsteady.yaml -p general:name:ALM_Unsteady_np_2 - - name: Run parallel regression tests within WindSE - shell: pwsh - run: | - conda activate test-environment - cd tests/9-Regression - mpirun -n 2 python -u ../../windse_driver/driver.py run ALM_Unsteady.yaml -p general:name:ALM_Unsteady_np_2 - # mpirun -n 2 python -u ../../windse_driver/driver.py run Floating_ALM.yaml -p general:name:Floating_ALM_np_2 + # - name: Run parallel regression tests within WindSE + # shell: pwsh + # run: | + # conda activate test-environment + # cd tests/9-Regression + # mpirun -n 2 python -u ../../windse_driver/driver.py run ALM_Unsteady.yaml -p general:name:ALM_Unsteady_np_2 + # # mpirun -n 2 python -u ../../windse_driver/driver.py run Floating_ALM.yaml -p general:name:Floating_ALM_np_2 - name: Run all regression tests within WindSE shell: pwsh diff --git a/.github/workflows/CI_every_PR_parallel.yml b/.github/workflows/CI_every_PR_parallel.yml new file mode 100644 index 00000000..78f510e5 --- /dev/null +++ b/.github/workflows/CI_every_PR_parallel.yml @@ -0,0 +1,50 @@ +name: CI_every_PR + +# We run CI on push commits on all branches +on: [pull_request] + +# A workflow run is made up of one or more jobs that can run sequentially or in parallel +jobs: + build: + name: CI_every_PR + runs-on: macos-latest + + steps: + - uses: actions/checkout@v2 + - uses: conda-incubator/setup-miniconda@v2 + # https://github.com/marketplace/actions/setup-miniconda + with: + miniconda-version: "latest" + channels: conda-forge + auto-update-conda: true + # python-version: 3.8 + # environment-file: environment.yml + + - name: Install test environment + shell: pwsh # putting in a shell command makes for compile linking problems later + # (if you use the shell here, cannot use 'compiler' package, but mpi only seems to work with it) + run: | + bash ./install.sh test-environment + conda activate test-environment + conda list + + # - name: Run parallel regression tests within WindSE + # shell: pwsh + # run: | + # conda activate test-environment + # pytest -sv --cov=windse tests/test_regression_parallel.py + + # - name: Run parallel regression tests within WindSE + # shell: pwsh + # run: | + # conda activate test-environment + # cd tests/9-Regression + # mpirun -n 2 windse run ALM_Unsteady.yaml -p general:name:ALM_Unsteady_np_2 + + - name: Run parallel regression tests within WindSE + shell: pwsh + run: | + conda activate test-environment + cd tests/9-Regression + mpirun -n 2 python -u ../../windse_driver/driver.py run ALM_Unsteady.yaml -p general:name:ALM_Unsteady_np_2 + # mpirun -n 2 python -u ../../windse_driver/driver.py run Floating_ALM.yaml -p general:name:Floating_ALM_np_2 diff --git a/tests/9-Regression/ALMx2_Unsteady.yaml b/tests/9-Regression/ALMx2_Unsteady.yaml index e8d5d6a0..d83b0f72 100644 --- a/tests/9-Regression/ALMx2_Unsteady.yaml +++ b/tests/9-Regression/ALMx2_Unsteady.yaml @@ -1,6 +1,6 @@ general: # name: "alm_original" - output: ["initial_guess","height","turbine_force","solution"] + output: ["initial_guess","height","turbine_force","solution","mesh"] # dolfin_adjoint: True debug_mode: True diff --git a/windse/ProblemManager.py b/windse/ProblemManager.py index 0fb20eca..8abc127b 100644 --- a/windse/ProblemManager.py +++ b/windse/ProblemManager.py @@ -79,12 +79,32 @@ def DebugOutput(self): else: e1 = Constant((1,0)); e2 = Constant((0,1)); - int_tf_x = assemble(inner(self.tf,e1)*dx)/self.dom.volume + if self.farm.use_local_tf_dx: + int_tf_x = 0 + int_tf_y = 0 + int_tf_z = 0 + for i in range(len(self.farm.tf_list)): + tf = self.farm.tf_list[i] + local_dx = self.farm.local_dx(i+1) + + int_tf_x += assemble(inner(tf,e1)*local_dx)/self.dom.volume + int_tf_y += assemble(inner(tf,e2)*local_dx)/self.dom.volume + if self.dom.dim == 3: + int_tf_z += assemble(inner(tf,e3)*local_dx)/self.dom.volume + else: + tf = 0 + for i in range(len(self.farm.tf_list)): + tf += self.farm.tf_list[i] + + int_tf_x = assemble(inner(tf,e1)*dx)/self.dom.volume + int_tf_y = assemble(inner(tf,e2)*dx)/self.dom.volume + if self.dom.dim == 3: + int_tf_z = assemble(inner(tf,e3)*dx)/self.dom.volume + self.tag_output("int_tf_x", int_tf_x) - int_tf_y = assemble(inner(self.tf,e2)*dx)/self.dom.volume self.tag_output("int_tf_y", int_tf_y) + if self.dom.dim == 3: - int_tf_z = assemble(inner(self.tf,e3)*dx)/self.dom.volume self.tag_output("int_tf_z", int_tf_z) if self.farm.turbine_type == 'line': @@ -110,25 +130,29 @@ def DebugOutput(self): self.tag_output("num_actuator_nodes", np.mean(num_actuator_nodes)) - def ComputeTurbineForce(self,u,inflow_angle,**kwargs): + def ComputeTurbineForceTerm(self,u,v,inflow_angle,**kwargs): tf_start = time.time() self.fprint("Calculating Turbine Force",special="header") - ### Compute the relative yaw angle ### + # Compute the relative yaw angle if inflow_angle is None: inflow_angle = self.dom.inflow_angle self.fprint('Computing turbine forces using %s' % (self.farm.turbine_type.upper())) - ### Create the turbine force function ### + # Create a list of turbine force function and domain of integration for each turbine if self.farm.turbine_type == "disabled" or self.farm.numturbs == 0: - tf = Function(self.fs.V) + + # if there are no turbine return an zero force term + tf_term = Function(self.fs.V)*dx else: - tf = self.farm.compute_turbine_force(u,inflow_angle,self.fs,**kwargs) + + # compute tf and dx for each turbine + tf_term = self.farm.compute_turbine_force(u,v,inflow_angle,self.fs,**kwargs) tf_stop = time.time() self.fprint("Turbine Force Calculated: {:1.2f} s".format(tf_stop-tf_start),special="footer") - return tf + return tf_term def ComputeTurbulenceModel(self, u): self.fprint(f"Using Turbulence Model: {self.turbulence_model}") @@ -294,7 +318,15 @@ def ComputeFunctional(self,inflow_angle): # mem0=memory_usage()[0] # mem_out, self.tf = memory_usage((self.ComputeTurbineForce,(self.u_k,inflow_angle),{}),max_usage=True,retval=True,max_iterations=1) # self.fprint("Memory Used: {:1.2f} MB".format(mem_out-mem0)) - self.tf = self.ComputeTurbineForce(self.u_k,inflow_angle) + tf_term = self.ComputeTurbineForceTerm(self.u_k,v,inflow_angle) + + # set_log_level(LogLevel.DEBUG) + # tic = time.time() + # assemble(turbine_force_term) + # toc = time.time() + # print(f"assemble time: {toc-tic} s") + # exit() + ### These constants will be moved into the params file ### f = Constant((0.0,)*self.dom.dim) @@ -321,7 +353,7 @@ def ComputeFunctional(self,inflow_angle): self.F += - inner(div(v),self.p_k)*dx self.F += - inner(div(self.u_k),q)*dx self.F += - inner(f,v)*dx - self.F += - inner(self.tf,v)*dx + self.F += - tf_term # Add body force to functional @@ -414,10 +446,16 @@ def ComputeFunctional(self,inflow_angle): self.nu_T=self.ComputeTurbulenceModel(self.u_k) ### Create the turbine force ### - self.tf = self.ComputeTurbineForce(self.u_k,inflow_angle) + tf_term = self.ComputeTurbineForceTerm(self.u_k,v,inflow_angle) ### Create the functional ### - self.F = inner(grad(self.u_k)*self.u_k, v)*dx + (nu+self.nu_T)*inner(grad(self.u_k), grad(v))*dx - inner(div(v),self.p_k)*dx - inner(div(self.u_k),q)*dx - inner(f,v)*dx - inner(self.tf,v)*dx + self.F = inner(grad(self.u_k)*self.u_k, v)*dx + self.F += (nu+self.nu_T)*inner(grad(self.u_k), grad(v))*dx + self.F += - inner(div(v),self.p_k)*dx + self.F += - inner(div(self.u_k),q)*dx + self.F += - inner(f,v)*dx + self.F += - tf_term + # Add body force to functional if abs(float(self.mbody_force)) >= 1e-14: @@ -512,7 +550,7 @@ def ComputeFunctional(self,inflow_angle): # self.tf = self.farm.TurbineForce(self.fs, self.dom.mesh, self.u_k2) # self.tf = Function(self.fs.V) - self.tf = self.ComputeTurbineForce(self.u_k,inflow_angle) + tf_term = self.ComputeTurbineForceTerm(self.u_k,v,inflow_angle) # self.u_k.assign(self.bd.bc_velocity) # self.u_k2.vector()[:] = 0.0 @@ -523,7 +561,7 @@ def ComputeFunctional(self,inflow_angle): # Solve for u_hat, a velocity estimate which doesn't include pressure gradient effects F1 = inner(dot(self.u_k, nabla_grad(u)), v)*dx \ + (nu+self.nu_T)*inner(grad(u), grad(v))*dx \ - - inner(self.tf, v)*dx + - tf_term # Add body force to functional F1 += inner(-self.mbody_force*self.bd.inflow_unit_vector,v)*dx @@ -543,7 +581,7 @@ def ComputeFunctional(self,inflow_angle): # Solve for u_star, a predicted velocity which includes the pressure gradient F3 = inner(dot(self.u_k, nabla_grad(u)), v)*dx \ + (nu+self.nu_T)*inner(grad(u), grad(v))*dx \ - - inner(self.tf, v)*dx \ + - tf_term \ + inner(grad(self.p_k), v)*dx \ + self.dt_1*inner(u - self.u_k, v)*dx @@ -666,7 +704,8 @@ def ComputeFunctional(self,inflow_angle): # self.tf = self.farm.TurbineForce(self.fs, self.dom.mesh, self.u_k2) # self.tf = Function(self.fs.V) - self.tf = self.ComputeTurbineForce(self.u_k,inflow_angle,simTime=0.0,simTime_prev=None, dt=self.dt) + tf_term = self.ComputeTurbineForceTerm(self.u_k,v,inflow_angle,simTime=0.0,simTime_prev=None, dt=self.dt) + self.u_k.assign(self.bd.bc_velocity) # Only the actuator lines point "upstream" against the flow @@ -699,7 +738,7 @@ def ComputeFunctional(self,inflow_angle): + inner(dot(U_AB, nabla_grad(U_CN)), v)*dx \ + (nu_c+self.nu_T)*inner(grad(U_CN), grad(v))*dx \ + inner(grad(self.p_k1), v)*dx \ - - dot(self.tf, v)*dx + - tf_term self.a1 = lhs(F1) self.L1 = rhs(F1) diff --git a/windse/SolverManager.py b/windse/SolverManager.py index 3a9a24c9..44f125da 100644 --- a/windse/SolverManager.py +++ b/windse/SolverManager.py @@ -878,8 +878,8 @@ def Solve(self): self.problem.bd.SaveInitialGuess(val=self.iter_val) if "height" in self.params.output and self.problem.dom.dim == 3: self.problem.bd.SaveHeight(val=self.iter_val) - # if "turbine_force" in self.params.output: - # self.problem.farm.SaveTurbineForce(val=self.iter_val) + if "turbine_force" in self.params.output: + self.problem.farm.save_functions(val=self.iter_val) self.fprint("Finished",special="footer") @@ -1110,9 +1110,9 @@ def save_adj(adj): # t1 = time.time() pr.enable() if 'dolfin' in self.problem.farm.turbines[0].type: - self.problem.ComputeTurbineForce(self.problem.u_k, self.problem.bd.inflow_angle, simTime=self.simTime, simTime_prev=self.simTime_prev, dt=self.problem.dt) + self.problem.ComputeTurbineForceTerm(self.problem.u_k, TestFunction(self.problem.fs.V), self.problem.bd.inflow_angle, simTime=self.simTime, simTime_prev=self.simTime_prev, dt=self.problem.dt) else: - # updated method from dev, is this necessary? + # updated method from dev, is this necessary? this may not work if there are self.farm.too_many_turbines self.problem.farm.update_turbine_force(self.problem.u_k, self.problem.bd.inflow_angle, self.problem.fs, simTime=self.simTime, simTime_prev=self.simTime_prev, dt=self.problem.dt) # new_tf = self.problem.ComputeTurbineForce(self.problem.u_k, self.problem.bd.inflow_angle, simTime=self.simTime, simTime_prev=self.simTime_prev, dt=self.problem.dt) # self.problem.tf.assign(new_tf) @@ -1316,10 +1316,10 @@ def SaveTimeSeries(self, simTime, simIter=None): self.problem.tf_save.vector()[:] = self.problem.tf_save.vector()[:] + fun.tf.vector()[:] else: if hasattr(self.problem,"tf_save"): - self.problem.tf_save = project(self.problem.tf, self.problem.fs.V, solver_type='cg') + self.problem.tf_save = project(sum(self.problem.farm.tf_list), self.problem.fs.V, solver_type='cg') else: self.problem.tf_save = Function(self.problem.fs.V) - self.problem.tf_save = project(self.problem.tf, self.problem.fs.V, solver_type='cg') + self.problem.tf_save = project(sum(self.problem.farm.tf_list), self.problem.fs.V, solver_type='cg') if self.first_save: diff --git a/windse/input_schema.yaml b/windse/input_schema.yaml index 9615376b..a43bc3aa 100644 --- a/windse/input_schema.yaml +++ b/windse/input_schema.yaml @@ -449,6 +449,18 @@ properties: - "null" description: "Seed to control/prescribe the randomness between runs." units: "dimensionless" + use_local_tf_dx: + default: True + type: + - "boolean" + description: "A flag to assemble the tf using a local dx measure for each turbine." + units: "dimensionless" + local_dx_scaling: + default: 2.5 + type: + - "number" + description: "Controls how big the local domain integration is in rotor diameters." + units: "dimensionless" # ================================================================ turbines: type: "object" diff --git a/windse/turbine_types/ActuatorDisk2D.py b/windse/turbine_types/ActuatorDisk2D.py index 64951574..1794daad 100644 --- a/windse/turbine_types/ActuatorDisk2D.py +++ b/windse/turbine_types/ActuatorDisk2D.py @@ -13,6 +13,8 @@ def __init__(self, i,x,y,dom,imported_params=None): def power(self, u, inflow_angle): x=SpatialCoordinate(self.dom.mesh) + print("here") + mx = self.mx my = self.my mz = self.mz diff --git a/windse/turbine_types/ActuatorLine.py b/windse/turbine_types/ActuatorLine.py index 2c48a1df..dcab809f 100644 --- a/windse/turbine_types/ActuatorLine.py +++ b/windse/turbine_types/ActuatorLine.py @@ -944,7 +944,7 @@ def power(self, u, inflow_angle): ys=self.my, zs=self.mz) - self.power_dolfin = assemble(1e-6*dot(-self.tf*self.angular_velocity, self.cyld_expr)*dx) + self.power_dolfin = 1e-6*dot(-self.tf*self.angular_velocity, self.cyld_expr) self.power_numpy = 1e-6*self.rotor_torque*self.angular_velocity # print("in turb.poweer()",self.power_dolfin, self.power_numpy) @@ -959,7 +959,7 @@ def prepare_saved_functions(self, func_list): [self.tf,"turbine_force"] ] else: - func_list[1][0] += self.tf + func_list[0][0] += self.tf return func_list diff --git a/windse/turbine_types/GenericTurbine.py b/windse/turbine_types/GenericTurbine.py index 380d88d5..c57fd8a1 100644 --- a/windse/turbine_types/GenericTurbine.py +++ b/windse/turbine_types/GenericTurbine.py @@ -5,7 +5,7 @@ from windse.blocks import blockify, BaseHeightBlock # import dolfin functions -from . import cos, sin, Constant +from . import cos, sin, Constant, MeshFunction, CompiledSubDomain, Measure # Other imports import warnings diff --git a/windse/wind_farm_types/GenericWindFarm.py b/windse/wind_farm_types/GenericWindFarm.py index ca69b4be..4b91b8af 100644 --- a/windse/wind_farm_types/GenericWindFarm.py +++ b/windse/wind_farm_types/GenericWindFarm.py @@ -2,7 +2,7 @@ from windse.turbine_types import turbine_dict import numpy as np import time, os -from . import MeshFunction, cells, project, FiniteElement, FunctionSpace, MixedElement, assemble, dx, parameters, Form, File +from . import MeshFunction, CompiledSubDomain, Measure, cells, project, inner, FiniteElement, FunctionSpace, MixedElement, assemble, dx, parameters, Form, File import matplotlib.pyplot as plt from pyadjoint.tape import stop_annotating from pyadjoint import AdjFloat @@ -28,6 +28,8 @@ def __init__(self, dom): # Load any global parameter self.farm_type = self.params["wind_farm"]["type"] + self.local_dx_scaling = self.params["wind_farm"]["local_dx_scaling"] + self.use_local_tf_dx = self.params["wind_farm"]["use_local_tf_dx"] self.turbine_type = self.params["turbines"]["type"] # Set any global flags @@ -52,7 +54,6 @@ def __init__(self, dom): ### Set the number of turbine when there are just too many and we need to project before solving self.too_many_turbines = 60 - # Setup the wind farm (eventually this will happen outside the init) self.fprint(f"Generating {self.name}",special="header") self.setup() @@ -169,97 +170,187 @@ def update_heights(self): for turb in self.turbines: turb.calculate_heights() - def compute_turbine_force(self,u,inflow_angle,fs,**kwargs): + def compute_turbine_force(self,u,v,inflow_angle,fs,**kwargs): """ Iterates over the turbines and adds up the turbine forces """ - # store the function space + # store the function space for later use self.fs = fs - # compute tf!!! - # TODO: make this more robust - if self.numturbs > self.too_many_turbines and "disk" in self.turbine_type: - # evenly distribute the number of turbine per project - n_groups = int((((self.numturbs-1)-((self.numturbs-1)%self.too_many_turbines))+self.too_many_turbines)/self.too_many_turbines) - split_value = round(self.numturbs/n_groups) + # create a list of turbine forces and integration domain + self.tf_list = [] + + # Create a cell function that will be used to store the subdomains + self.turbine_subdomains = MeshFunction("size_t", self.dom.mesh, self.dom.mesh.topology().dim()) + self.turbine_subdomains.set_all(0) + + # Iterate over each turbine to build its individual subdomain + for turb in self.turbines: + self.tf_list.append(turb.turbine_force(u,inflow_angle,fs,**kwargs)) + + # the subdomain is a sphere centered at the hub of the turbine with a height of scaling*RD and diameter of scaling*RD + + if self.dom.dim == 3: + subdomain_marker = CompiledSubDomain("(x[0]-x0)*(x[0]-x0)+(x[1]-y0)*(x[1]-y0)+(x[2]-HH)*(x[2]-HH)<=r*r",r=self.local_dx_scaling*turb.RD/2.0, HH=turb.z, x0=turb.x, y0=turb.y) + else: + subdomain_marker = CompiledSubDomain("(x[0]-x0)*(x[0]-x0)+(x[1]-y0)*(x[1]-y0)<=r*r",r=self.local_dx_scaling*turb.RD/2.0, x0=turb.x, y0=turb.y) + + + # add the subdomain to the marker + subdomain_marker.mark(self.turbine_subdomains,turb.index+1) + + + # create a measure for the turbines + self.local_dx = Measure('dx', subdomain_data=self.turbine_subdomains) + + # check if the number of turbines might trigger a python recursion error. + if self.numturbs > self.too_many_turbines: + final_tf_list, final_local_dx = self.compress_turbine_function(u) + else: + final_tf_list = self.tf_list + final_local_dx = self.local_dx + + # Combine the force and dx to create a single tf term + if self.use_local_tf_dx: + tf_term = 0 + for i in range(len(final_tf_list)): + tf_term += inner(final_tf_list[i],v)*final_local_dx(i+1) + + # use the full domain of integration if local is not requested + else: + tf_term = inner(sum(final_tf_list),v)*dx + + return tf_term + + def compress_turbine_function(self, u): + """ + used to assemble the turbine force in chunks to avoid recursion errors + """ + + # evenly distribute the number of turbine per project + n_groups = int((((self.numturbs-1)-((self.numturbs-1)%self.too_many_turbines))+self.too_many_turbines)/self.too_many_turbines) + split_value = round(self.numturbs/n_groups) + print(f"Oh dear, there are just too many turbines (>={self.too_many_turbines}). I'm just going to collect them into {n_groups} groups of around {split_value} turbines each and then project.") + + # init loop values + combined_tf_list = [] + combined_dx_tf_list = [] + counter = 0 + group = 1 + + # Create a cell function that will be used to store the subdomains + temp_subdomain = MeshFunction("size_t", self.dom.mesh, self.dom.mesh.topology().dim()) + temp_subdomain.set_all(0) - print(f"Oh dear, there are just too many turbines (>={self.too_many_turbines}). I'm just going to collect them into {n_groups} groups of around {split_value} turbines each and then project.") + # disks need to be treated specially due to how tf is convoluted with u + if "disk" in self.turbine_type: - # init loop values - tf1_proj = [] - tf2_proj = [] - tf3_proj = [] + # init combined tf tf1 = 0 tf2 = 0 tf3 = 0 - counter = 0 - group = 1 + + # sum and project tf parts for turb in self.turbines: - turb.turbine_force(u,inflow_angle,fs,**kwargs) + + # get the turbine parts tf1 += turb.tf1 tf2 += turb.tf2 tf3 += turb.tf3 + + # the subdomain is a sphere centered at the hub of the turbine with a height of scaling*RD and diameter of scaling*RD + if self.dom.dim == 3: + subdomain_marker = CompiledSubDomain("(x[0]-x0)*(x[0]-x0)+(x[1]-y0)*(x[1]-y0)+(x[2]-HH)*(x[2]-HH)<=r*r",r=self.local_dx_scaling*turb.RD/2.0, HH=turb.z, x0=turb.x, y0=turb.y) + else: + subdomain_marker = CompiledSubDomain("(x[0]-x0)*(x[0]-x0)+(x[1]-y0)*(x[1]-y0)<=r*r",r=self.local_dx_scaling*turb.RD/2.0, x0=turb.x, y0=turb.y) + + # add the subdomain to the marker + subdomain_marker.mark(temp_subdomain,group) + + # once we have collected the maximum number of turbines, combine them an reset the sums if counter >= split_value: self.fprint(f"Projecting tf1 for group {group}") - tf1_proj.append(project(tf1,fs.V,solver_type='gmres',preconditioner_type="hypre_amg")) + proj_tf1 = project(tf1,fs.V,solver_type='gmres',preconditioner_type="hypre_amg") self.fprint(f"Projecting tf2 for group {group}") - tf2_proj.append(project(tf2,fs.V,solver_type='gmres',preconditioner_type="hypre_amg")) + proj_tf2 = project(tf2,fs.V,solver_type='gmres',preconditioner_type="hypre_amg") self.fprint(f"Projecting tf3 for group {group}") - tf3_proj.append(project(tf3,fs.V,solver_type='gmres',preconditioner_type="hypre_amg")) + proj_tf3 = project(tf3,fs.V,solver_type='gmres',preconditioner_type="hypre_amg") + combined_tf_list.append(proj_tf1*u[0]**2+proj_tf2*u[1]**2+proj_tf3*u[0]*u[1]) + + # reset loop values tf1 = 0 tf2 = 0 tf3 = 0 counter = 0 group += 1 + temp_subdomain = MeshFunction("size_t", self.dom.mesh, self.dom.mesh.topology().dim()) + temp_subdomain.set_all(0) + else: counter += 1 + # project the final turbine group self.fprint(f"Projecting tf1 for group {group}") - tf1_proj.append(project(tf1,fs.V,solver_type='gmres',preconditioner_type="hypre_amg")) + proj_tf1 = project(tf1,fs.V,solver_type='gmres',preconditioner_type="hypre_amg") self.fprint(f"Projecting tf2 for group {group}") - tf2_proj.append(project(tf2,fs.V,solver_type='gmres',preconditioner_type="hypre_amg")) + proj_tf2 = project(tf2,fs.V,solver_type='gmres',preconditioner_type="hypre_amg") self.fprint(f"Projecting tf3 for group {group}") - tf3_proj.append(project(tf3,fs.V,solver_type='gmres',preconditioner_type="hypre_amg")) - - - # do one more sum and project - tf1 = sum(tf1_proj) - tf2 = sum(tf2_proj) - tf3 = sum(tf3_proj) + proj_tf3 = project(tf3,fs.V,solver_type='gmres',preconditioner_type="hypre_amg") + combined_tf_list.append(proj_tf1*u[0]**2+proj_tf2*u[1]**2+proj_tf3*u[0]*u[1]) - # self.fprint(f"Final project tf1") - # tf1 = project(sum(tf1_proj),fs.V) - # self.fprint(f"Final project tf2") - # tf2 = project(sum(tf2_proj),fs.V) - # self.fprint(f"Final project tf3") - # tf3 = project(sum(tf3_proj),fs.V) + # create a measure for the final turbine group + final_local_dx.append(Measure('dx', subdomain_data=temp_subdomain)) - # create the final turbine force - return tf1*u[0]**2+tf2*u[1]**2+tf3*u[0]*u[1] + # combine turbine for every other style of turbine else: - # do the normal loop + + # init combined tf tf = 0 + + # sum and project tf parts for turb in self.turbines: - tf += turb.turbine_force(u,inflow_angle,fs,**kwargs) - return tf - - # T = 0 - # dt = 0.5 - # max_its = 20 - # tf_file = File("output/float_demo/pitch_tf.pvd") - # for i in range(max_its): - # print(f"Generating TF at time: {T} s") - # tf = 0 - # for turb in self.turbines: - # tf += turb.turbine_force(u,inflow_angle,fs,simTime=T) - # tf = project(tf,fs.V,solver_type='cg',preconditioner_type="hypre_amg") - # tf.rename("tf","tf") - # tf_file << (tf,T) - # T += dt - # exit() - # return tf + # get the turbine parts + tf += turb.tf + + # the subdomain is a sphere centered at the hub of the turbine with a height of scaling*RD and diameter of scaling*RD + if self.dom.dim == 3: + subdomain_marker = CompiledSubDomain("(x[0]-x0)*(x[0]-x0)+(x[1]-y0)*(x[1]-y0)+(x[2]-HH)*(x[2]-HH)<=r*r",r=self.local_dx_scaling*turb.RD/2.0, HH=turb.z, x0=turb.x, y0=turb.y) + else: + subdomain_marker = CompiledSubDomain("(x[0]-x0)*(x[0]-x0)+(x[1]-y0)*(x[1]-y0)<=r*r",r=self.local_dx_scaling*turb.RD/2.0, x0=turb.x, y0=turb.y) + + # add the subdomain to the marker + subdomain_marker.mark(temp_subdomain,group) + + # once we have collected the maximum number of turbines, combine them an reset the sums + if counter >= split_value: + self.fprint(f"Projecting tf for group {group}") + proj_tf = project(tf,fs.V,solver_type='gmres',preconditioner_type="hypre_amg") + combined_tf_list.append(proj_tf) + + # create a measure for this turbine group + combined_dx_tf_list.append(Measure('dx', subdomain_data=temp_subdomain)) + + # reset loop values + tf = 0 + counter = 0 + group += 1 + temp_subdomain = MeshFunction("size_t", self.dom.mesh, self.dom.mesh.topology().dim()) + temp_subdomain.set_all(0) + else: + counter += 1 + + # project the final turbine group + self.fprint(f"Projecting tf for group {group}") + proj_tf = project(tf,fs.V,solver_type='gmres',preconditioner_type="hypre_amg") + combined_tf_list.append(proj_tf) + + # create a measure for the final turbine group + final_local_dx.append(Measure('dx', subdomain_data=temp_subdomain)) + + # create the final turbine force + return combined_tf_list, final_local_dx def update_controls(self): """ @@ -283,11 +374,18 @@ def compute_power(self, u, inflow_angle): ### Build the integrand val = 0 for turb in self.turbines: - val += turb.power(u, inflow_angle) + temp = turb.power(u, inflow_angle) + if not isinstance(val,(float,AdjFloat)): + if self.use_local_tf_dx: + val += temp*self.local_dx(turb.index+1) + else: + val += temp*dx + else: + val += temp ### Assemble if needed if not isinstance(val,(float,AdjFloat)): - J = assemble(val*dx) + J = assemble(val) else: J = val @@ -308,7 +406,10 @@ def save_power(self, u, inflow_angle, iter_val = 0.0, simTime = 0.0): ### Assemble if needed if not isinstance(val,(float,AdjFloat)): - J = assemble(val*dx) + if self.use_local_tf_dx: + J = assemble(val*self.local_dx(turb.index+1)) + else: + J = assemble(val*dx) else: J = val @@ -501,11 +602,11 @@ def save_functions(self,val=0): # gather functions func_list = [] - if self.numturbs>=self.too_many_turbines: - print(f"Oh dear, there are just too many turbines (>={self.too_many_turbines}). Saving turbine function is disabled.") - else: - for turb in self.turbines: - func_list = turb.prepare_saved_functions(func_list) + # if self.numturbs>=self.too_many_turbines: + # print(f"Oh dear, there are just too many turbines (>={self.too_many_turbines}). Saving turbine function is disabled.") + # else: + for turb in self.turbines: + func_list = turb.prepare_saved_functions(func_list) # prepare to store file pointer if needed if self.func_first_save: @@ -535,6 +636,12 @@ def save_functions(self,val=0): else: self.params.Save(func, func_name,subfolder="functions/",val=val,file=self.func_files[i]) + # save the farm level turbine subdomain function + if self.func_first_save: + self.subdomain_file = self.params.Save(self.turbine_subdomains,"turbine_subdomains",subfolder="mesh/",val=val,filetype="pvd") + else: + self.params.Save(self.turbine_subdomains,"turbine_subdomains",subfolder="mesh/",val=val,file=self.subdomain_file,filetype="pvd") + # Flip the flag! self.func_first_save = False diff --git a/windse_driver/driver.py b/windse_driver/driver.py index 5fb59f1b..c492bd6c 100644 --- a/windse_driver/driver.py +++ b/windse_driver/driver.py @@ -110,6 +110,9 @@ def run_action(params_loc=None): params.comm.Barrier() + + dolfin.list_timings(dolfin.TimingClear.clear, [dolfin.TimingType.wall]) + return runtime def mesh_action(params_loc=None): From 7e9ffcc7915e1c24d52a1eb2f7d6fce968fb605b Mon Sep 17 00:00:00 2001 From: jefalon Date: Fri, 1 Mar 2024 16:15:45 -0700 Subject: [PATCH 02/15] changed name of new parallel CI file --- .github/workflows/CI_every_PR_parallel.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI_every_PR_parallel.yml b/.github/workflows/CI_every_PR_parallel.yml index 78f510e5..4231da53 100644 --- a/.github/workflows/CI_every_PR_parallel.yml +++ b/.github/workflows/CI_every_PR_parallel.yml @@ -1,4 +1,4 @@ -name: CI_every_PR +name: CI_every_PR_parallel # We run CI on push commits on all branches on: [pull_request] From 192e4a2b223ec9462c18d400eef6ebf5e8a2f7ee Mon Sep 17 00:00:00 2001 From: jefalon Date: Tue, 12 Mar 2024 11:28:50 -0600 Subject: [PATCH 03/15] added local dx test problem --- .github/workflows/CI_every_PR_parallel.yml | 2 +- .../9-Regression/3D_Box_Steady_Local_dx.yaml | 83 ++++++++++++ .../3D_Box_Steady_Local_dx_truth.yaml | 121 ++++++++++++++++++ windse/ProblemManager.py | 41 ++++-- windse/input_schema.yaml | 2 +- windse/turbine_types/ActuatorDisk2D.py | 2 - windse/wind_farm_types/GenericWindFarm.py | 20 ++- 7 files changed, 243 insertions(+), 28 deletions(-) create mode 100644 tests/9-Regression/3D_Box_Steady_Local_dx.yaml create mode 100644 tests/9-Regression/Truth_Data/3D_Box_Steady_Local_dx_truth.yaml diff --git a/.github/workflows/CI_every_PR_parallel.yml b/.github/workflows/CI_every_PR_parallel.yml index 4231da53..9d844e61 100644 --- a/.github/workflows/CI_every_PR_parallel.yml +++ b/.github/workflows/CI_every_PR_parallel.yml @@ -6,7 +6,7 @@ on: [pull_request] # A workflow run is made up of one or more jobs that can run sequentially or in parallel jobs: build: - name: CI_every_PR + name: CI_every_PR_parallel runs-on: macos-latest steps: diff --git a/tests/9-Regression/3D_Box_Steady_Local_dx.yaml b/tests/9-Regression/3D_Box_Steady_Local_dx.yaml new file mode 100644 index 00000000..8ef4d737 --- /dev/null +++ b/tests/9-Regression/3D_Box_Steady_Local_dx.yaml @@ -0,0 +1,83 @@ +# General options +general: + # name: 3D_Wind_Farm # Name of the output folder + output: ["mesh","initial_guess","height","turbine_force","solution"] + dolfin_adjoint: True + debug_mode: True + +# Wind Farm Parameters: +wind_farm: + + ########################## Grid Wind Farm ######################### + type: grid # | + grid_rows: 2 # Number of rows | - + grid_cols: 4 # Number of columns | - + ex_x: [-600, 600] # x-extent of the farm | m + ex_y: [-300, 300] # y-extent of the farm | m + use_local_tf_dx: True + local_dx_scaling: 2.5 + ################################################################### + +turbines: + type: disk # Use actuator Disks | + HH: 90 # Hub Height | m + RD: 126.0 # Turbine Diameter | m + thickness: 20.0 # Effective Thickness | m + yaw: 0.349 # Yaw | rads + axial: 0.33 # Axial Induction | - + +domain: + + ########################### Box Domain ############################ + type: box # | + x_range: [-1200, 1200] # x-range of the domain | m + y_range: [-600, 600] # y-range of the domain | m + z_range: [0.04, 640] # z-range of the domain | m + nx: 16 # Number of x-nodes | - + ny: 8 # Number of y-nodes | - + nz: 8 # Number of z-nodes | - + ################################################################### + +refine: + warp_type: split + warp_percent: 0.7 # percent of cells moved | - + warp_height: 250 # move cell below this value | m + refine_custom: + 1: + type: box + x_range: [-1000,1200] + y_range: [-400, 400] + z_range: [0, 200] + + turbine_num: 1 # number of turbine refinements| - + turbine_factor: 1.25 # turbine radius multiplier | - + +function_space: + type: linear + +boundary_conditions: + vel_profile: log + HH_vel: 8.0 + k: 0.4 + # inflow_angle: 1.13 + ######## Uncomment to test out custom BCs for BoxDomain ########## + boundary_types: + inflow: ["west"] + no_stress: ["east"] + free_slip: ["top","north","south"] + no_slip: ["bottom"] + ################################################################## + +problem: + type: stabilized + viscosity: 5 + lmax: 50 + +solver: + type: steady + save_power: true + +optimization: + control_types: [layout,yaw] + layout_bounds: [[-720, 720],[-450, 450]] + gradient: True diff --git a/tests/9-Regression/Truth_Data/3D_Box_Steady_Local_dx_truth.yaml b/tests/9-Regression/Truth_Data/3D_Box_Steady_Local_dx_truth.yaml new file mode 100644 index 00000000..97a02375 --- /dev/null +++ b/tests/9-Regression/Truth_Data/3D_Box_Steady_Local_dx_truth.yaml @@ -0,0 +1,121 @@ +GenericWindFarm: + min_x: -537.0 + max_x: 537.0 + avg_x: 0.0 + min_y: -237.0 + max_y: 237.0 + avg_y: 0.0 + min_z: 90.0 + max_z: 90.0 + avg_z: 90.0 + min_yaw: 0.349 + max_yaw: 0.349 + avg_yaw: 0.349 +DomainManager: + dim: 3 + x_min: -1200.0 + x_max: 1200.0 + y_min: -600.0 + y_max: 600.0 + z_min: 0.04 + z_max: 640.0 + volume: 1843084799.9992278 + ID_east: 1 + SA_east: 767951.9999999994 + ID_north: 2 + SA_north: 1535904.0000000016 + ID_west: 3 + SA_west: 767952.0000000009 + ID_south: 4 + SA_south: 1535904.0000000016 + ID_bottom: 5 + SA_bottom: 2880000.0000000005 + ID_top: 6 + SA_top: 2880000.000000001 +FunctionSpaceManager: + velocity_dofs: 66945 + pressure_dofs: 22315 + total_dofs: 89260 +BoundaryManager: + min_x: 3.4975322130308117 + max_x: 9.14464606512413 + avg_x: 7.815257995353568 + min_y: 0.0 + max_y: 0.0 + avg_y: 0.0 + min_z: 0.0 + max_z: 0.0 + avg_z: 0.0 + min_p: 0.0 + max_p: 0.0 + avg_p: 0.0 + min_initial_values: 0.0 + max_initial_values: 9.14464606512413 + avg_initial_values: 1.953814498838392 + num_bc: 11 +ProblemManager: + int_nu_T: 2.331420005441958 + int_tf_x: -0.002864684105259432 + int_tf_y: -0.001042446119382597 + int_tf_z: 0.0 +SolverManager: + min_x_vel: 0.0 + max_x_vel: 10.0110040797921 + avg_x_vel: 8.328121562200876 + min_y_vel: -1.1608289338114017 + max_y_vel: 0.8388368562010192 + avg_y_vel: 0.01134324259391201 + min_z_vel: -0.7269189939175986 + max_z_vel: 0.8452329732622681 + avg_z_vel: 0.01134324259391201 +OptimizationManager: + n_controls: 24 + obj_value0: 6.382825298587173 + val0_x_0: -537.0 + val0_y_0: -237.0 + val0_x_1: -179.0 + val0_y_1: -237.0 + val0_x_2: 179.0 + val0_y_2: -237.0 + val0_x_3: 537.0 + val0_y_3: -237.0 + val0_x_4: -537.0 + val0_y_4: 237.0 + val0_x_5: -179.0 + val0_y_5: 237.0 + val0_x_6: 179.0 + val0_y_6: 237.0 + val0_x_7: 537.0 + val0_y_7: 237.0 + val0_yaw_0: 0.349 + val0_yaw_1: 0.349 + val0_yaw_2: 0.349 + val0_yaw_3: 0.349 + val0_yaw_4: 0.349 + val0_yaw_5: 0.349 + val0_yaw_6: 0.349 + val0_yaw_7: 0.349 + grad_x_0: -0.002944577519889316 + grad_y_0: -0.003879559290036203 + grad_x_1: 0.0002572885785810967 + grad_y_1: -0.000771343276552843 + grad_x_2: 0.00016125200033123392 + grad_y_2: 0.0009408947252452374 + grad_x_3: 0.0006326386029324466 + grad_y_3: 0.001092719658152443 + grad_x_4: -0.001087667520347604 + grad_y_4: -0.00041629066649610715 + grad_x_5: -0.0003968572825411281 + grad_y_5: -0.0006407825878008103 + grad_x_6: -0.0004920257466021231 + grad_y_6: 0.0009641726599597918 + grad_x_7: 0.0013265487291185697 + grad_y_7: 0.0022613991905987233 + grad_yaw_0: -0.568027762149109 + grad_yaw_1: -0.3264574252231357 + grad_yaw_2: -0.3422881185400642 + grad_yaw_3: -0.45540462187448805 + grad_yaw_4: -0.6715042516148744 + grad_yaw_5: -0.36733320180969986 + grad_yaw_6: -0.3528333077534691 + grad_yaw_7: -0.4720037582351555 diff --git a/windse/ProblemManager.py b/windse/ProblemManager.py index 8abc127b..e8c585c4 100644 --- a/windse/ProblemManager.py +++ b/windse/ProblemManager.py @@ -416,6 +416,9 @@ class TaylorHoodProblem(GenericProblem): def __init__(self,domain,windfarm,function_space,boundary_conditions): super(TaylorHoodProblem, self).__init__(domain,windfarm,function_space,boundary_conditions) + + self.first_loop = True + ### Create Functional ### self.ComputeFunctional(self.dom.inflow_angle) self.DebugOutput() @@ -426,41 +429,53 @@ def ComputeFunctional(self,inflow_angle): ### These constants will be moved into the params file ### f = Constant((0.0,)*self.dom.dim) vonKarman=0.41 - eps=Constant(0.01) + eps=Constant(0.000001) nu = self.viscosity self.fprint("Viscosity: {:1.2e}".format(float(self.viscosity))) self.fprint("Max Mixing Length: {:1.2e}".format(float(self.lmax))) + ### Create the test/trial/functions ### - self.up_k = Function(self.fs.W) - self.u_k,self.p_k = split(self.up_k) - v,q = TestFunctions(self.fs.W) + self.v,self.q = TestFunctions(self.fs.W) ### Set the initial guess ### ### (this will become a separate function.) + # if self.first_loop: + # self.up_k = Function(self.fs.W) + # self.up_k.assign(self.bd.u0) + # self.first_loop = False + # else: + # temp_up = self.up_k.copy() + self.up_k = Function(self.fs.W) self.up_k.assign(self.bd.u0) + self.u_k,self.p_k = split(self.up_k) ### Calculate nu_T self.nu_T=self.ComputeTurbulenceModel(self.u_k) ### Create the turbine force ### - tf_term = self.ComputeTurbineForceTerm(self.u_k,v,inflow_angle) + tf_term = self.ComputeTurbineForceTerm(self.u_k,self.v,inflow_angle) ### Create the functional ### - self.F = inner(grad(self.u_k)*self.u_k, v)*dx - self.F += (nu+self.nu_T)*inner(grad(self.u_k), grad(v))*dx - self.F += - inner(div(v),self.p_k)*dx - self.F += - inner(div(self.u_k),q)*dx - self.F += - inner(f,v)*dx + self.F = inner(grad(self.u_k)*self.u_k, self.v)*dx + self.F += (nu+self.nu_T)*inner(grad(self.u_k), grad(self.v))*dx + self.F += - inner(div(self.v),self.p_k)*dx + self.F += - inner(div(self.u_k),self.q)*dx + # self.F += - inner(f,v)*dx self.F += - tf_term + # ### Add in the Stabilizing term ### + # if abs(float(eps)) >= 1e-14: + # self.fprint("Using Stabilization Term") + # stab = - eps*inner(grad(q), grad(self.p_k))*dx - eps*inner(grad(q), dot(grad(self.u_k), self.u_k))*dx + # self.F += stab # Add body force to functional if abs(float(self.mbody_force)) >= 1e-14: self.fprint("Using Body Force") - self.F += inner(-self.mbody_force*self.bd.inflow_unit_vector,v)*dx + self.F += inner(-self.mbody_force*self.bd.inflow_unit_vector,self.v)*dx if self.use_25d_model: if self.dom.dim == 3: @@ -471,9 +486,9 @@ def ComputeFunctional(self,inflow_angle): dvdy = Dx(self.u_k[1], 1) if inflow_angle is None: - term25 = dvdy*q*dx + term25 = dvdy*self.q*dx else: - term25 = (abs(sin(inflow_angle))*dudx*q + abs(cos(inflow_angle))*dvdy*q)*dx + term25 = (abs(sin(inflow_angle))*dudx*self.q + abs(cos(inflow_angle))*dvdy*self.q)*dx self.F -= term25 diff --git a/windse/input_schema.yaml b/windse/input_schema.yaml index a43bc3aa..18291dd0 100644 --- a/windse/input_schema.yaml +++ b/windse/input_schema.yaml @@ -450,7 +450,7 @@ properties: description: "Seed to control/prescribe the randomness between runs." units: "dimensionless" use_local_tf_dx: - default: True + default: False type: - "boolean" description: "A flag to assemble the tf using a local dx measure for each turbine." diff --git a/windse/turbine_types/ActuatorDisk2D.py b/windse/turbine_types/ActuatorDisk2D.py index 1794daad..64951574 100644 --- a/windse/turbine_types/ActuatorDisk2D.py +++ b/windse/turbine_types/ActuatorDisk2D.py @@ -13,8 +13,6 @@ def __init__(self, i,x,y,dom,imported_params=None): def power(self, u, inflow_angle): x=SpatialCoordinate(self.dom.mesh) - print("here") - mx = self.mx my = self.my mz = self.mz diff --git a/windse/wind_farm_types/GenericWindFarm.py b/windse/wind_farm_types/GenericWindFarm.py index 4b91b8af..af396858 100644 --- a/windse/wind_farm_types/GenericWindFarm.py +++ b/windse/wind_farm_types/GenericWindFarm.py @@ -205,7 +205,7 @@ def compute_turbine_force(self,u,v,inflow_angle,fs,**kwargs): # check if the number of turbines might trigger a python recursion error. if self.numturbs > self.too_many_turbines: - final_tf_list, final_local_dx = self.compress_turbine_function(u) + final_tf_list, final_local_dx = self.compress_turbine_function(u,fs) else: final_tf_list = self.tf_list final_local_dx = self.local_dx @@ -222,7 +222,7 @@ def compute_turbine_force(self,u,v,inflow_angle,fs,**kwargs): return tf_term - def compress_turbine_function(self, u): + def compress_turbine_function(self, u, fs): """ used to assemble the turbine force in chunks to avoid recursion errors """ @@ -234,7 +234,7 @@ def compress_turbine_function(self, u): # init loop values combined_tf_list = [] - combined_dx_tf_list = [] + final_local_dx = [] counter = 0 group = 1 @@ -284,8 +284,6 @@ def compress_turbine_function(self, u): tf3 = 0 counter = 0 group += 1 - temp_subdomain = MeshFunction("size_t", self.dom.mesh, self.dom.mesh.topology().dim()) - temp_subdomain.set_all(0) else: counter += 1 @@ -300,7 +298,7 @@ def compress_turbine_function(self, u): combined_tf_list.append(proj_tf1*u[0]**2+proj_tf2*u[1]**2+proj_tf3*u[0]*u[1]) # create a measure for the final turbine group - final_local_dx.append(Measure('dx', subdomain_data=temp_subdomain)) + final_local_dx = Measure('dx', subdomain_data=temp_subdomain) # combine turbine for every other style of turbine else: @@ -328,9 +326,6 @@ def compress_turbine_function(self, u): self.fprint(f"Projecting tf for group {group}") proj_tf = project(tf,fs.V,solver_type='gmres',preconditioner_type="hypre_amg") combined_tf_list.append(proj_tf) - - # create a measure for this turbine group - combined_dx_tf_list.append(Measure('dx', subdomain_data=temp_subdomain)) # reset loop values tf = 0 @@ -347,7 +342,7 @@ def compress_turbine_function(self, u): combined_tf_list.append(proj_tf) # create a measure for the final turbine group - final_local_dx.append(Measure('dx', subdomain_data=temp_subdomain)) + final_local_dx = Measure('dx', subdomain_data=temp_subdomain) # create the final turbine force return combined_tf_list, final_local_dx @@ -401,7 +396,6 @@ def save_power(self, u, inflow_angle, iter_val = 0.0, simTime = 0.0): J_list[1]=simTime with stop_annotating(): for i,turb in enumerate(self.turbines): - val = turb.power(u,inflow_angle) ### Assemble if needed @@ -600,6 +594,10 @@ def save_functions(self,val=0): Note: this function is way over engineered! """ + if self.numturbs >= self.too_many_turbines: + print(f"Oh dear, there are just too many turbines so saving the turbine force is disabled") + return + # gather functions func_list = [] # if self.numturbs>=self.too_many_turbines: From 464faba840c65573d7796576a8a1fc717f11db78 Mon Sep 17 00:00:00 2001 From: jefalon Date: Tue, 12 Mar 2024 12:44:09 -0600 Subject: [PATCH 04/15] fixed typo in power calc --- windse/wind_farm_types/GenericWindFarm.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/windse/wind_farm_types/GenericWindFarm.py b/windse/wind_farm_types/GenericWindFarm.py index af396858..2edb735b 100644 --- a/windse/wind_farm_types/GenericWindFarm.py +++ b/windse/wind_farm_types/GenericWindFarm.py @@ -370,7 +370,7 @@ def compute_power(self, u, inflow_angle): val = 0 for turb in self.turbines: temp = turb.power(u, inflow_angle) - if not isinstance(val,(float,AdjFloat)): + if not isinstance(temp,(float,AdjFloat)): if self.use_local_tf_dx: val += temp*self.local_dx(turb.index+1) else: From 1690ac351ca708497488915a818121f1c24dcc93 Mon Sep 17 00:00:00 2001 From: jefalon Date: Wed, 8 May 2024 12:43:57 -0600 Subject: [PATCH 05/15] refactor wind inflow angle to be consistent with meteorlogical convention --- demo/documented/Driver_Example/parameters.rst | 2 +- .../Yaml_Examples/0-wind_farm_2D.yaml | 2 +- .../1-wind_farm_2D_layout_opt.yaml | 2 +- .../Yaml_Examples/2-wind_farm_3D.yaml | 4 +-- .../Yaml_Examples/3-multiangle_solve.yaml | 4 +-- .../4-multiangle_optimization.yaml | 2 +- .../Yaml_Examples/5-yaw_optimization.yaml | 2 +- .../Yaml_Examples/7-wind_farm_2D_OM.yaml | 2 +- .../Yaml_Examples/8-spacing-and-shear.yaml | 2 +- .../9-save_objective_output.yaml | 2 +- .../Input_Data/iea_3p4/wind_farm.csv | 2 +- .../Input_Data/nrel_5mw/wind_farm.csv | 6 ++-- doc/source/param_tips.rst | 7 +++++ doc/source/params_old.rst | 2 +- .../2D_Rectangle_Steady_gmsh.yaml | 2 +- tests/9-Regression/3D_Box_Steady.yaml | 2 +- .../9-Regression/3D_Box_Steady_Local_dx.yaml | 2 +- tests/9-Regression/3D_Box_Steady_gmsh.yaml | 2 +- tests/9-Regression/3D_Skew_Steady.yaml | 4 +-- .../Input_Data/ALM_Unsteady/iea_rwt.csv | 2 +- .../ALMx2_Unsteady/nrel_5mw/wind_farm.csv | 4 +-- .../offset_farms/offset_turbines_4.csv | 2 +- .../yawed_farms/yawed_turbines_0.csv | 4 +-- .../yawed_farms/yawed_turbines_1.csv | 2 +- .../yawed_farms/yawed_turbines_2.csv | 2 +- .../yawed_farms/yawed_turbines_3.csv | 2 +- .../yawed_farms/yawed_turbines_4.csv | 2 +- .../yawed_farms/yawed_turbines_5.csv | 2 +- .../yawed_farms/yawed_turbines_6.csv | 2 +- .../yawed_farms/yawed_turbines_7.csv | 2 +- .../Yaw_Optimization/Farm_Layout.csv | 4 +-- windse/BoundaryManager.py | 1 + windse/DomainManager.py | 18 ++++++++---- windse/OptimizationManager.py | 15 +++++++++- windse/PostprocessingManager.py | 2 +- windse/SolverManager.py | 28 +++++++++++-------- windse/helper_functions.py | 24 ++++++++++++++++ windse/input_schema.yaml | 8 +++--- windse/turbine_types/ActuatorDisk.py | 9 +++--- windse/turbine_types/ActuatorDisk2D.py | 2 +- windse/turbine_types/ActuatorDiskExpr.py | 5 ++-- windse/turbine_types/ActuatorHybridDisk.py | 5 ++-- windse/turbine_types/GenericTurbine.py | 7 +++++ windse/wind_farm_types/GenericWindFarm.py | 2 +- 44 files changed, 138 insertions(+), 71 deletions(-) diff --git a/demo/documented/Driver_Example/parameters.rst b/demo/documented/Driver_Example/parameters.rst index ddfe0a7a..ba76aad8 100755 --- a/demo/documented/Driver_Example/parameters.rst +++ b/demo/documented/Driver_Example/parameters.rst @@ -76,7 +76,7 @@ the rotor diameter to ensure all turbines including the rotors are located within the extents. The rest of the parameters determine the physical properties of the turbines: - * ``yaw``: The yaw of the turbines where 0 is perpendicular to an East to West inflow. + * ``yaw``: The yaw of the turbines where 0 is perpendicular to inflow. * ``axial``: The axial induction * ``HH``: The hub height relative to the ground * ``RD``: The rotor diameter diff --git a/demo/documented/Yaml_Examples/0-wind_farm_2D.yaml b/demo/documented/Yaml_Examples/0-wind_farm_2D.yaml index af7f84c6..ba9fc3f0 100755 --- a/demo/documented/Yaml_Examples/0-wind_farm_2D.yaml +++ b/demo/documented/Yaml_Examples/0-wind_farm_2D.yaml @@ -43,7 +43,7 @@ turbines: # These two are left off because they are supplied in the wind_farm.csv file # HH: 90 # Hub Height | m - # yaw: 0.0 # Yaw | rads + # yaw: 0.0 # Yaw | degs ################################################################### # Domain Parameters: Uncomment a set to change domain shape diff --git a/demo/documented/Yaml_Examples/1-wind_farm_2D_layout_opt.yaml b/demo/documented/Yaml_Examples/1-wind_farm_2D_layout_opt.yaml index 3b52b4f4..6ff22c0d 100644 --- a/demo/documented/Yaml_Examples/1-wind_farm_2D_layout_opt.yaml +++ b/demo/documented/Yaml_Examples/1-wind_farm_2D_layout_opt.yaml @@ -40,7 +40,7 @@ turbines: HH: 90 # Hub Height | m RD: 126.0 # Turbine Diameter | m thickness: 12.0 # Effective Thickness | m - yaw: 0.0 # Yaw | rads + yaw: 0.0 # Yaw | degs axial: 0.25 # Axial Induction | - # Domain Parameters: Uncomment a set to change domain shape diff --git a/demo/documented/Yaml_Examples/2-wind_farm_3D.yaml b/demo/documented/Yaml_Examples/2-wind_farm_3D.yaml index 71b8bee2..da1d65a8 100644 --- a/demo/documented/Yaml_Examples/2-wind_farm_3D.yaml +++ b/demo/documented/Yaml_Examples/2-wind_farm_3D.yaml @@ -23,7 +23,7 @@ turbines: # These two are left off because they are supplied in the wind_farm.csv file # HH: 90 # Hub Height | m - # yaw: 0.0 # Yaw | rads + # yaw: 0.0 # Yaw | degs ################################################################### # Domain Parameters: Uncomment a set to change domain shape @@ -78,7 +78,7 @@ boundary_conditions: vel_profile: log HH_vel: 8.0 k: 0.4 - inflow_angle: 1.13 + inflow_angle: 285.0 ######### Uncomment to test out custom BCs for BoxDomain ########## # boundary_types: # inflow: ["west"] diff --git a/demo/documented/Yaml_Examples/3-multiangle_solve.yaml b/demo/documented/Yaml_Examples/3-multiangle_solve.yaml index e525a4d0..4f6bee16 100755 --- a/demo/documented/Yaml_Examples/3-multiangle_solve.yaml +++ b/demo/documented/Yaml_Examples/3-multiangle_solve.yaml @@ -23,7 +23,7 @@ turbines: HH: 90 # Hub Height | m RD: 126.0 # Turbine Diameter | m thickness: 20.0 # Effective Thickness | m - yaw: 0.0 # Yaw | rads + yaw: 0.0 # Yaw | degs axial: 0.33 # Axial Induction | - # Domain Parameters: Uncomment a set to change domain shape @@ -64,7 +64,7 @@ boundary_conditions: vel_profile: log HH_vel: 8.0 k: 0.4 - inflow_angle: [1.34,5.04,3] + inflow_angle: [45,115,3] boundary_types: inflow: ["west"] no_stress: ["east"] diff --git a/demo/documented/Yaml_Examples/4-multiangle_optimization.yaml b/demo/documented/Yaml_Examples/4-multiangle_optimization.yaml index c0edc1df..2e8d4295 100644 --- a/demo/documented/Yaml_Examples/4-multiangle_optimization.yaml +++ b/demo/documented/Yaml_Examples/4-multiangle_optimization.yaml @@ -21,7 +21,7 @@ turbines: HH: 80 # Hub Height | m RD: 80.0 # Turbine Diameter | m thickness: 20 # Effective Thickness | m - yaw: 0.0 # Yaw | rads + yaw: 0.0 # Yaw | degs axial: 0.33 # Axial Induction | - # Domain Constants for a Box Domain diff --git a/demo/documented/Yaml_Examples/5-yaw_optimization.yaml b/demo/documented/Yaml_Examples/5-yaw_optimization.yaml index 7bd07a9f..6f8f1165 100755 --- a/demo/documented/Yaml_Examples/5-yaw_optimization.yaml +++ b/demo/documented/Yaml_Examples/5-yaw_optimization.yaml @@ -24,7 +24,7 @@ turbines: # These two are left off because they are supplied in the wind_farm.csv file # HH: 90 # Hub Height | m - # yaw: 0.0 # Yaw | rads + # yaw: 0.0 # Yaw | degs ################################################################### diff --git a/demo/documented/Yaml_Examples/7-wind_farm_2D_OM.yaml b/demo/documented/Yaml_Examples/7-wind_farm_2D_OM.yaml index 8f51ac8e..aac6167b 100644 --- a/demo/documented/Yaml_Examples/7-wind_farm_2D_OM.yaml +++ b/demo/documented/Yaml_Examples/7-wind_farm_2D_OM.yaml @@ -40,7 +40,7 @@ turbines: HH: 90 # Hub Height | m RD: 126.0 # Turbine Diameter | m thickness: 12.0 # Effective Thickness | m - yaw: 0.0 # Yaw | rads + yaw: 0.0 # Yaw | degs axial: 0.25 # Axial Induction | - # Domain Parameters: Uncomment a set to change domain shape diff --git a/demo/documented/Yaml_Examples/8-spacing-and-shear.yaml b/demo/documented/Yaml_Examples/8-spacing-and-shear.yaml index 25a0c5c7..522914f0 100644 --- a/demo/documented/Yaml_Examples/8-spacing-and-shear.yaml +++ b/demo/documented/Yaml_Examples/8-spacing-and-shear.yaml @@ -26,7 +26,7 @@ turbines: HH: 90 # Hub Height | m RD: 126.0 # Turbine Diameter | m thickness: 12.0 # Effective Thickness | m - yaw: 0.0 # Yaw | rads + yaw: 0.0 # Yaw | degs axial: 0.25 # Axial Induction | - # Domain Parameters: Uncomment a set to change domain shape diff --git a/demo/documented/Yaml_Examples/9-save_objective_output.yaml b/demo/documented/Yaml_Examples/9-save_objective_output.yaml index 21880c2d..9dca76dd 100644 --- a/demo/documented/Yaml_Examples/9-save_objective_output.yaml +++ b/demo/documented/Yaml_Examples/9-save_objective_output.yaml @@ -15,7 +15,7 @@ turbines: HH: 90 RD: 126.0 thickness: 20.0 - yaw: 0.349 + yaw: 20.0 axial: 0.33 domain: diff --git a/demo/documented/Yaml_Examples/Input_Data/iea_3p4/wind_farm.csv b/demo/documented/Yaml_Examples/Input_Data/iea_3p4/wind_farm.csv index 80460a2a..cd3576cb 100644 --- a/demo/documented/Yaml_Examples/Input_Data/iea_3p4/wind_farm.csv +++ b/demo/documented/Yaml_Examples/Input_Data/iea_3p4/wind_farm.csv @@ -1,2 +1,2 @@ x, y, HH, Yaw, Diameter, Thickness, Axial_Induction -0, 0, 110.0, -0.349, 130.0, 13.0, .33 +0, 0, 110.0, -20.0, 130.0, 13.0, .33 diff --git a/demo/documented/Yaml_Examples/Input_Data/nrel_5mw/wind_farm.csv b/demo/documented/Yaml_Examples/Input_Data/nrel_5mw/wind_farm.csv index da2d6586..094ca3b0 100644 --- a/demo/documented/Yaml_Examples/Input_Data/nrel_5mw/wind_farm.csv +++ b/demo/documented/Yaml_Examples/Input_Data/nrel_5mw/wind_farm.csv @@ -1,4 +1,4 @@ x, y, HH, yaw, Diameter, Thickness, Axial_Induction - -882, 0, 90.0, 0.262, 126.0, 12.0, 0.25 - 0, 0, 90.0, 0.000, 126.0, 12.0, 0.25 - 882, 0, 90.0, -0.262, 126.0, 12.0, 0.25 \ No newline at end of file + -882, 0, 90.0, 15.0, 126.0, 12.0, 0.25 + 0, 0, 90.0, 0.0, 126.0, 12.0, 0.25 + 882, 0, 90.0, -15.0, 126.0, 12.0, 0.25 \ No newline at end of file diff --git a/doc/source/param_tips.rst b/doc/source/param_tips.rst index ea7a240c..6f416c93 100644 --- a/doc/source/param_tips.rst +++ b/doc/source/param_tips.rst @@ -203,6 +203,13 @@ The syntax for each refinement type is:: +.. _inflow_angle: + +Defining Inflow Angle +--------------------- + +The inflow angle is defined in the meteorological sense with units of degrees. This means an inflow angle of 0 will result in wind blowing from the north to the south. Increasing the inflow angle will rotate the wind clockwise (i.e. 90 degrees corresponds to wind blowing from the east and moving west). The default value of 270 degrees, which corresponds to wind blowing from the west to the east. + .. _custom_boundaries: Customizing Boundary Conditions diff --git a/doc/source/params_old.rst b/doc/source/params_old.rst index 8e2077d9..f072d976 100644 --- a/doc/source/params_old.rst +++ b/doc/source/params_old.rst @@ -359,7 +359,7 @@ This section will define all the parameters for the wind farm:: | ``RD`` | The rotor diameter | all | None | m | | | | | | | +------------------------+-----------------------------------------------+--------------------+----------+-------------+ -| ``yaw`` | | Determines the yaw of all turbines. Yaw is | all | None | rad | +| ``yaw`` | | Determines the yaw of all turbines. Yaw is | all | None | degs | | | | relative to the wind inflow direction | | | | +------------------------+-----------------------------------------------+--------------------+----------+-------------+ | ``thickness`` | The effective thickness of the rotor disk | | "disk" or disk | None | m | diff --git a/tests/9-Regression/2D_Rectangle_Steady_gmsh.yaml b/tests/9-Regression/2D_Rectangle_Steady_gmsh.yaml index 3652e39d..a8500727 100644 --- a/tests/9-Regression/2D_Rectangle_Steady_gmsh.yaml +++ b/tests/9-Regression/2D_Rectangle_Steady_gmsh.yaml @@ -21,7 +21,7 @@ turbines: HH: 90 # Hub Height | m RD: 126.0 # Turbine Diameter | m thickness: 20.0 # Effective Thickness | m - yaw: 0.349 # Yaw | rads + yaw: 19.996 # Yaw | degs axial: 0.33 # Axial Induction | - domain: diff --git a/tests/9-Regression/3D_Box_Steady.yaml b/tests/9-Regression/3D_Box_Steady.yaml index ab9370ad..a7cd9bf9 100644 --- a/tests/9-Regression/3D_Box_Steady.yaml +++ b/tests/9-Regression/3D_Box_Steady.yaml @@ -21,7 +21,7 @@ turbines: HH: 90 # Hub Height | m RD: 126.0 # Turbine Diameter | m thickness: 20.0 # Effective Thickness | m - yaw: 0.349 # Yaw | rads + yaw: 19.996 # Yaw | degs axial: 0.33 # Axial Induction | - domain: diff --git a/tests/9-Regression/3D_Box_Steady_Local_dx.yaml b/tests/9-Regression/3D_Box_Steady_Local_dx.yaml index 8ef4d737..3d69c89f 100644 --- a/tests/9-Regression/3D_Box_Steady_Local_dx.yaml +++ b/tests/9-Regression/3D_Box_Steady_Local_dx.yaml @@ -23,7 +23,7 @@ turbines: HH: 90 # Hub Height | m RD: 126.0 # Turbine Diameter | m thickness: 20.0 # Effective Thickness | m - yaw: 0.349 # Yaw | rads + yaw: 19.996 # Yaw | degs axial: 0.33 # Axial Induction | - domain: diff --git a/tests/9-Regression/3D_Box_Steady_gmsh.yaml b/tests/9-Regression/3D_Box_Steady_gmsh.yaml index 210c39f4..112fa919 100644 --- a/tests/9-Regression/3D_Box_Steady_gmsh.yaml +++ b/tests/9-Regression/3D_Box_Steady_gmsh.yaml @@ -21,7 +21,7 @@ turbines: HH: 90 # Hub Height | m RD: 126.0 # Turbine Diameter | m thickness: 20.0 # Effective Thickness | m - yaw: 0.349 # Yaw | rads + yaw: 19.996 # Yaw | degs axial: 0.33 # Axial Induction | - domain: diff --git a/tests/9-Regression/3D_Skew_Steady.yaml b/tests/9-Regression/3D_Skew_Steady.yaml index ab8d40b1..5f539d7d 100644 --- a/tests/9-Regression/3D_Skew_Steady.yaml +++ b/tests/9-Regression/3D_Skew_Steady.yaml @@ -21,7 +21,7 @@ turbines: HH: 90 # Hub Height | m RD: 126.0 # Turbine Diameter | m thickness: 10.5 # Effective Thickness | m - yaw: 0.0 # Yaw | rads + yaw: 0.0 # Yaw | degs axial: 0.33 # Axial Induction | - @@ -64,7 +64,7 @@ boundary_conditions: vel_profile: log HH_vel: 8.0 k: 0.4 - inflow_angle: 1.13 + inflow_angle: 205.256 problem: type: stabilized diff --git a/tests/9-Regression/Input_Data/ALM_Unsteady/iea_rwt.csv b/tests/9-Regression/Input_Data/ALM_Unsteady/iea_rwt.csv index 1c3c6f61..8e807a32 100644 --- a/tests/9-Regression/Input_Data/ALM_Unsteady/iea_rwt.csv +++ b/tests/9-Regression/Input_Data/ALM_Unsteady/iea_rwt.csv @@ -1,2 +1,2 @@ x, y, HH, Yaw, Diameter -0.0, 0.0, 110.0, -0.349, 130.0 +0.0, 0.0, 110.0, -19.996, 130.0 diff --git a/tests/9-Regression/Input_Data/ALMx2_Unsteady/nrel_5mw/wind_farm.csv b/tests/9-Regression/Input_Data/ALMx2_Unsteady/nrel_5mw/wind_farm.csv index da2d6586..0e807f62 100644 --- a/tests/9-Regression/Input_Data/ALMx2_Unsteady/nrel_5mw/wind_farm.csv +++ b/tests/9-Regression/Input_Data/ALMx2_Unsteady/nrel_5mw/wind_farm.csv @@ -1,4 +1,4 @@ x, y, HH, yaw, Diameter, Thickness, Axial_Induction - -882, 0, 90.0, 0.262, 126.0, 12.0, 0.25 + -882, 0, 90.0, 15.011, 126.0, 12.0, 0.25 0, 0, 90.0, 0.000, 126.0, 12.0, 0.25 - 882, 0, 90.0, -0.262, 126.0, 12.0, 0.25 \ No newline at end of file + 882, 0, 90.0, -15.011, 126.0, 12.0, 0.25 \ No newline at end of file diff --git a/tests/9-Regression/Input_Data/ALMx2_Unsteady/offset_farms/offset_turbines_4.csv b/tests/9-Regression/Input_Data/ALMx2_Unsteady/offset_farms/offset_turbines_4.csv index 748d5295..5ac75a63 100644 --- a/tests/9-Regression/Input_Data/ALMx2_Unsteady/offset_farms/offset_turbines_4.csv +++ b/tests/9-Regression/Input_Data/ALMx2_Unsteady/offset_farms/offset_turbines_4.csv @@ -1,3 +1,3 @@ x, y, HH, yaw 0, 0, 0.0, 0.000 - 630, 72, 0.0, 0.000 + 630, 72, 0.0, 0.000 diff --git a/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_0.csv b/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_0.csv index ad9c9ff5..51bb9222 100644 --- a/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_0.csv +++ b/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_0.csv @@ -1,3 +1,3 @@ x, y, HH, yaw - 0, 0, 0.0, 0.000 - 630, 0, 0.0, 0.000 + 0, 0, 0.0, 270.000 + 630, 0, 0.0, 270.000 diff --git a/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_1.csv b/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_1.csv index 4a9a6af2..bba20fbc 100644 --- a/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_1.csv +++ b/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_1.csv @@ -1,3 +1,3 @@ x, y, HH, yaw - 0, 0, 0.0, -0.08726646 + 0, 0, 0.0, -5.0 630, 0, 0.0, 0.000 diff --git a/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_2.csv b/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_2.csv index 290605e7..561a8de7 100644 --- a/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_2.csv +++ b/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_2.csv @@ -1,3 +1,3 @@ x, y, HH, yaw - 0, 0, 0.0, -0.17453293 + 0, 0, 0.0, -10.0 630, 0, 0.0, 0.000 diff --git a/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_3.csv b/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_3.csv index ab2e21c7..1e5d042b 100644 --- a/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_3.csv +++ b/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_3.csv @@ -1,3 +1,3 @@ x, y, HH, yaw - 0, 0, 0.0, -0.26179939 + 0, 0, 0.0, -15.0 630, 0, 0.0, 0.000 diff --git a/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_4.csv b/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_4.csv index 447a8c3e..bf940edb 100644 --- a/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_4.csv +++ b/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_4.csv @@ -1,3 +1,3 @@ x, y, HH, yaw - 0, 0, 0.0, -0.34906585 + 0, 0, 0.0, -20.0 630, 0, 0.0, 0.000 diff --git a/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_5.csv b/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_5.csv index 71a0e145..6f55e831 100644 --- a/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_5.csv +++ b/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_5.csv @@ -1,3 +1,3 @@ x, y, HH, yaw - 0, 0, 0.0, -0.43633231 + 0, 0, 0.0, -25.0 630, 0, 0.0, 0.000 diff --git a/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_6.csv b/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_6.csv index 0b701e3c..8247c3ec 100644 --- a/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_6.csv +++ b/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_6.csv @@ -1,3 +1,3 @@ x, y, HH, yaw - 0, 0, 0.0, -0.52359878 + 0, 0, 0.0, -30.0 630, 0, 0.0, 0.000 diff --git a/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_7.csv b/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_7.csv index 4edfe92a..3ef7c24e 100644 --- a/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_7.csv +++ b/tests/9-Regression/Input_Data/Hybrid_Disk/yawed_farms/yawed_turbines_7.csv @@ -1,3 +1,3 @@ x, y, HH, yaw - 0, 0, 0.0, -0.61086524 + 0, 0, 0.0, -35.0 630, 0, 0.0, 0.000 diff --git a/tests/9-Regression/Input_Data/Yaw_Optimization/Farm_Layout.csv b/tests/9-Regression/Input_Data/Yaw_Optimization/Farm_Layout.csv index 6b167c55..65600251 100644 --- a/tests/9-Regression/Input_Data/Yaw_Optimization/Farm_Layout.csv +++ b/tests/9-Regression/Input_Data/Yaw_Optimization/Farm_Layout.csv @@ -1,3 +1,3 @@ x, y, HH, yaw, RD, thickness, axial -0.0, -325.0, 110.0, 0.5236, 130.0, 13.0, 0.33 -0.0, 325.0, 110.0, -0.5236, 130.0, 13.0, 0.33 \ No newline at end of file +0.0, -325.0, 110.0, 30.0, 130.0, 13.0, 0.33 +0.0, 325.0, 110.0, -30.0, 130.0, 13.0, 0.33 \ No newline at end of file diff --git a/windse/BoundaryManager.py b/windse/BoundaryManager.py index 085995bc..7256fa49 100644 --- a/windse/BoundaryManager.py +++ b/windse/BoundaryManager.py @@ -44,6 +44,7 @@ def __init__(self,dom,fs,farm): ### Update attributes based on params file ### for key, value in self.params["boundary_conditions"].items(): setattr(self,key,value) + del self.inflow_angle # remove inflow angle from bc because domain owns it. self.extra_kwarg = {} if self.params.dolfin_adjoint: diff --git a/windse/DomainManager.py b/windse/DomainManager.py index e14e538a..416fe435 100644 --- a/windse/DomainManager.py +++ b/windse/DomainManager.py @@ -30,6 +30,9 @@ ### Import the cumulative parameters ### from windse import windse_parameters + ### Import some helper functions + from windse.helper_functions import meteor_to_math + ### Check if we need dolfin_adjoint ### if windse_parameters.dolfin_adjoint: from dolfin_adjoint import * @@ -187,11 +190,16 @@ def __init__(self): self.ground_reference = self.ground_reference*self.xscale ### Get the initial wind direction ### - self.inflow_angle = self.params["boundary_conditions"]["inflow_angle"] - if self.inflow_angle is None: - self.inflow_angle = 0.0 - elif isinstance(self.inflow_angle,list): - self.inflow_angle = self.inflow_angle[0] + self.raw_inflow_angle = self.params["boundary_conditions"]["inflow_angle"] + if self.raw_inflow_angle is None: + inflow_angle = 270.0 # default inflow angle is from west to east + elif isinstance(self.raw_inflow_angle,list): + inflow_angle = self.raw_inflow_angle[0] + else: + inflow_angle = self.raw_inflow_angle + print(inflow_angle) + self.inflow_angle = meteor_to_math(inflow_angle) + print(self.inflow_angle) self.initial_inflow_angle = self.inflow_angle ### Add a space to store periodic boundary conditions ### diff --git a/windse/OptimizationManager.py b/windse/OptimizationManager.py index f8d9aa20..8a5a8842 100644 --- a/windse/OptimizationManager.py +++ b/windse/OptimizationManager.py @@ -592,6 +592,11 @@ def Gradient(self): # d_format = np.float64(d) # self.params.comm.Gather(d_format, d_global, root=0) # d_sum = np.sum(d_global) + + if "yaw" in self.names[i]: + ctl_val = np.degrees(ctl_val) + d = d*np.pi/180 + d_out = '%12s: %12.5e, %22.15e' % (self.names[i], ctl_val, d) # print('Rank %d, %s' % (self.params.rank, d_out)) @@ -630,6 +635,8 @@ def ListControls(self,m): self.fprint("Next Control Values",special="header") for i, val in enumerate(m): + if "yaw" in self.names[i]: + val = np.degrees(float(val)) self.fprint(self.names[i] +": " +repr(float(val))) self.fprint("",special="footer") @@ -666,7 +673,13 @@ def SaveControls(self,m): else: m_old.append(float(getattr(self.farm.turbines[index],control_name))) # print(m_new) - # print(type(m_new)) + # print(type(m_new)) + + for i in range(len(m_new)): + if "yaw" in self.names[i]: + m_new[i] = np.degrees(float(m_new[i])) + m_old[i] = np.degrees(float(m_old[i])) + if self.problem.params.rank == 0: if self.iteration == 0: diff --git a/windse/PostprocessingManager.py b/windse/PostprocessingManager.py index 8b4bac19..7fa11388 100644 --- a/windse/PostprocessingManager.py +++ b/windse/PostprocessingManager.py @@ -81,7 +81,7 @@ def write_to_floris(data_to_write, solver): yaws = solver.problem.farm.get_yaw_angles() - floris_dict["farm"]["yaw"] = yaws + floris_dict["farm"]["yaw"] = np.degrees(yaws) # Generate the full path to the FLORIS file path_to_floris_file = os.path.join(solver.params.folder, diff --git a/windse/SolverManager.py b/windse/SolverManager.py index 44f125da..9e5f6916 100644 --- a/windse/SolverManager.py +++ b/windse/SolverManager.py @@ -29,6 +29,9 @@ ### Import the cumulative parameters ### from windse import windse_parameters + ### Import some helper functions + from windse.helper_functions import meteor_to_math, math_to_meteor + ### Check if we need dolfin_adjoint ### if windse_parameters.dolfin_adjoint: from dolfin_adjoint import * @@ -1110,11 +1113,11 @@ def save_adj(adj): # t1 = time.time() pr.enable() if 'dolfin' in self.problem.farm.turbines[0].type: - self.problem.ComputeTurbineForceTerm(self.problem.u_k, TestFunction(self.problem.fs.V), self.problem.bd.inflow_angle, simTime=self.simTime, simTime_prev=self.simTime_prev, dt=self.problem.dt) + self.problem.ComputeTurbineForceTerm(self.problem.u_k, TestFunction(self.problem.fs.V), self.problem.dom.inflow_angle, simTime=self.simTime, simTime_prev=self.simTime_prev, dt=self.problem.dt) else: # updated method from dev, is this necessary? this may not work if there are self.farm.too_many_turbines - self.problem.farm.update_turbine_force(self.problem.u_k, self.problem.bd.inflow_angle, self.problem.fs, simTime=self.simTime, simTime_prev=self.simTime_prev, dt=self.problem.dt) - # new_tf = self.problem.ComputeTurbineForce(self.problem.u_k, self.problem.bd.inflow_angle, simTime=self.simTime, simTime_prev=self.simTime_prev, dt=self.problem.dt) + self.problem.farm.update_turbine_force(self.problem.u_k, self.problem.dom.inflow_angle, self.problem.fs, simTime=self.simTime, simTime_prev=self.simTime_prev, dt=self.problem.dt) + # new_tf = self.problem.ComputeTurbineForce(self.problem.u_k, self.problem.dom.inflow_angle, simTime=self.simTime, simTime_prev=self.simTime_prev, dt=self.problem.dt) # self.problem.tf.assign(new_tf) pr.disable() @@ -1128,7 +1131,7 @@ def save_adj(adj): # # self.problem.alm_power_sum += self.problem.alm_power*self.problem.dt # # self.problem.alm_power_average = self.problem.alm_power_sum/self.simTime - # self.problem.alm_power_dolfin = self.problem.farm.compute_power(self.problem.u_k, self.problem.bd.inflow_angle) + # self.problem.alm_power_dolfin = self.problem.farm.compute_power(self.problem.u_k, self.problem.dom.inflow_angle) # output_str = 'Rotor Power (dolfin, solver): %s MW' % (self.problem.alm_power_dolfin/1.0e6) # self.fprint(output_str) @@ -2041,27 +2044,28 @@ def __init__(self,problem): if self.params["domain"]["type"] in ["imported"]: raise ValueError("Cannot use a Multi-Angle Solver with an "+self.params["domain"]["type"]+" domain.") self.orignal_solve = super(MultiAngleSolver, self).Solve - if self.problem.bd.inflow_angle is None: + if self.problem.dom.raw_inflow_angle is None: self.wind_range = [0, 2.0*np.pi,self.num_wind_angles] - elif isinstance(self.problem.bd.inflow_angle,list): - if len(self.problem.bd.inflow_angle)==3: - self.wind_range = self.problem.bd.inflow_angle + elif isinstance(self.problem.dom.raw_inflow_angle,list): + if len(self.problem.dom.raw_inflow_angle)==3: + self.wind_range = self.problem.dom.raw_inflow_angle else: - self.wind_range = [self.problem.bd.inflow_angle[0],self.problem.bd.inflow_angle[1],self.num_wind_angles] + self.wind_range = [self.problem.dom.raw_inflow_angle[0],self.problem.dom.raw_inflow_angle[1],self.num_wind_angles] else: - self.wind_range = [self.problem.bd.inflow_angle,self.problem.bd.inflow_angle+2.0*np.pi,self.num_wind_angles] + self.wind_range = [self.problem.dom.raw_inflow_angle,self.problem.dom.raw_inflow_angle+2.0*np.pi,self.num_wind_angles] self.angles = np.linspace(*self.wind_range,endpoint=self.endpoint) + self.angles = meteor_to_math(self.angles) # self.angles += self.angle_offset def Solve(self): for i, theta in enumerate(self.angles): self.fprint("Performing Solve {:d} of {:d}".format(i+1,len(self.angles)),special="header") - self.fprint("Wind Angle: "+repr(theta)) + self.fprint("Wind Angle: "+repr(math_to_meteor(theta))) if i > 0 or not near(theta,self.problem.dom.inflow_angle): self.problem.dom.inflow_angle = theta self.ChangeWindAngle(theta) - self.iter_val = theta + self.iter_val = math_to_meteor(theta) self.orignal_solve() self.fprint("Finished Solve {:d} of {:d}".format(i+1,len(self.angles)),special="footer") diff --git a/windse/helper_functions.py b/windse/helper_functions.py index 72859803..29e30717 100644 --- a/windse/helper_functions.py +++ b/windse/helper_functions.py @@ -9,6 +9,30 @@ from windse.blocks import blockify, MpiEvalBlock, UflEvalBlock from pyadjoint.overloaded_type import create_overloaded_object +def meteor_to_math(angle): + ''' + this function converts an angle from the meteorological convention to the math convention. + 0 -> 270 + 90 -> 180 + 180 -> 90 + 270 -> 0 + Additionally the input is in degrees and the output is in radians. + + ''' + return np.radians(270.0-angle) + +def math_to_meteor(angle): + ''' + this function converts an angle from the math convention to the meteorological convention. + 0 -> 270 + 90 -> 180 + 180 -> 90 + 270 -> 0 + Additionally the input is in radians and the output is in degrees. + + ''' + return 270.0-np.degrees(angle) + def ufl_eval(form, print_statement=None): ''' This function converts complex ufl forms to floats diff --git a/windse/input_schema.yaml b/windse/input_schema.yaml index d89c1c9f..67a6411c 100644 --- a/windse/input_schema.yaml +++ b/windse/input_schema.yaml @@ -503,7 +503,7 @@ properties: - "number" - "null" description: "Yaw of the turbine relative to inflow angle." - units: "radians" + units: "degrees" axial: default: Null type: @@ -769,12 +769,12 @@ properties: - "null" description: "Location for the turbsim inflow data." inflow_angle: - default: 0.0 + default: 270.0 type: - "number" - "array" - description: "Angle of the inflow velocity." - units: "radians" + description: "Angle of the inflow velocity. See :ref:`inflow_angle`" + units: "degrees" boundary_names: type: "object" description: "Used for renaming the boundries. See :ref:`custom_boundaries`" diff --git a/windse/turbine_types/ActuatorDisk.py b/windse/turbine_types/ActuatorDisk.py index d397119b..fbdf81f1 100644 --- a/windse/turbine_types/ActuatorDisk.py +++ b/windse/turbine_types/ActuatorDisk.py @@ -36,11 +36,11 @@ def create_controls(self): # The control list is very important for extracting information from the optimzation step. self.controls_list = ["x","y","yaw","axial"] - def build_actuator_disk(self,inflow_angle): + def build_actuator_disk(self): ### Alias useful values ### x0 = [self.mx,self.my,self.mz] - yaw = self.myaw+inflow_angle + yaw = self.myaw+self.inflow_angle W = self.thickness*1.0 R = self.RD/2.0 ma = self.maxial @@ -124,10 +124,11 @@ def turbine_force(self,u,inflow_angle,fs): ### compute the space kernal and radial force self.fs = fs - self.actuator_disk = self.build_actuator_disk(inflow_angle) + self.inflow_angle.assign(inflow_angle) + self.actuator_disk = self.build_actuator_disk() ### Expand the dot product - yaw = self.myaw+inflow_angle + yaw = self.myaw+self.inflow_angle self.tf1 = self.actuator_disk * cos(yaw)**2 self.tf2 = self.actuator_disk * sin(yaw)**2 self.tf3 = self.actuator_disk * 2.0 * cos(yaw) * sin(yaw) diff --git a/windse/turbine_types/ActuatorDisk2D.py b/windse/turbine_types/ActuatorDisk2D.py index 64951574..8e58da50 100644 --- a/windse/turbine_types/ActuatorDisk2D.py +++ b/windse/turbine_types/ActuatorDisk2D.py @@ -20,7 +20,7 @@ def power(self, u, inflow_angle): W = self.thickness*1.0 R = self.RD/2.0 ma = self.maxial - yaw = self.myaw+inflow_angle + yaw = self.myaw+self.inflow_angle A = np.pi*R**2.0 C_tprime = 4*ma/(1-ma) C_pprime = 0.45/(1-ma)**3 diff --git a/windse/turbine_types/ActuatorDiskExpr.py b/windse/turbine_types/ActuatorDiskExpr.py index 8bdceb69..6a253084 100644 --- a/windse/turbine_types/ActuatorDiskExpr.py +++ b/windse/turbine_types/ActuatorDiskExpr.py @@ -117,10 +117,11 @@ def turbine_force(self,u,inflow_angle,fs): """ ### compute the space kernal and radial force - self.actuator_disk = self.build_actuator_disk(inflow_angle) + self.inflow_angle.assign(inflow_angle) + self.actuator_disk = self.build_actuator_disk(self.inflow_angle) ### Expand the dot product - yaw = self.myaw+inflow_angle + yaw = self.myaw+self.inflow_angle tf1 = self.actuator_disk * cos(yaw)**2 tf2 = self.actuator_disk * sin(yaw)**2 tf3 = self.actuator_disk * 2.0 * cos(yaw) * sin(yaw) diff --git a/windse/turbine_types/ActuatorHybridDisk.py b/windse/turbine_types/ActuatorHybridDisk.py index 1b54928f..08908a59 100644 --- a/windse/turbine_types/ActuatorHybridDisk.py +++ b/windse/turbine_types/ActuatorHybridDisk.py @@ -98,7 +98,7 @@ def build_actuator_disk(self,inflow_angle): ### Alias useful values ### x0 = [self.mx,self.my,self.mz] - yaw = self.myaw+inflow_angle + yaw = self.myaw+self.inflow_angle W = self.thickness*1.0 R = self.RD/2.0 ma = self.maxial @@ -214,6 +214,7 @@ def turbine_force(self,u,inflow_angle,fs): ### Store function space ### self.fs = fs + self.inflow_angle.assign(inflow_angle) ### Initialize some blade parameters ### self.init_blade_properties() @@ -222,7 +223,7 @@ def turbine_force(self,u,inflow_angle,fs): self.actuator_disk = self.build_actuator_disk(inflow_angle) ### Expand the dot product - yaw = self.myaw+inflow_angle + yaw = self.myaw+self.inflow_angle self.tf1 = self.actuator_disk * cos(yaw)**2 self.tf2 = self.actuator_disk * sin(yaw)**2 self.tf3 = self.actuator_disk * 2.0 * cos(yaw) * sin(yaw) diff --git a/windse/turbine_types/GenericTurbine.py b/windse/turbine_types/GenericTurbine.py index c57fd8a1..70a74ea4 100644 --- a/windse/turbine_types/GenericTurbine.py +++ b/windse/turbine_types/GenericTurbine.py @@ -35,6 +35,7 @@ def __init__(self, i,x,y,dom,imported_params): self.x = x self.y = y self.imported_params = imported_params + self.inflow_angle = Constant(0.0) # blockify custom functions so dolfin adjoint can track them if self.params.performing_opt_calc: @@ -49,6 +50,7 @@ def setup(self): """ self.load_parameters() self.update_parameters_from_file() + self.convert_yaw() self.compute_parameters() self.create_controls() self.calculate_heights() @@ -87,6 +89,11 @@ def update_parameters_from_file(self): # update the parameter setattr(self,column_header,val) + def convert_yaw(self): + if hasattr(self,"yaw"): + self.yaw = np.radians(self.yaw) + + def compute_parameters(self): """ This function will compute any additional parameters diff --git a/windse/wind_farm_types/GenericWindFarm.py b/windse/wind_farm_types/GenericWindFarm.py index 2edb735b..505ec210 100644 --- a/windse/wind_farm_types/GenericWindFarm.py +++ b/windse/wind_farm_types/GenericWindFarm.py @@ -663,7 +663,7 @@ def save_wind_farm(self,val=None,filename="wind_farm"): ### Get quantities ### x, y, z = self.get_hub_locations().T HH = z-self.get_ground_heights() - yaw = self.get_yaw_angles() + yaw = np.degrees(self.get_yaw_angles()) RD = self.get_rotor_diameters() ### Might not need these anymore From 39146018fb2584d55f575ae08284a54916307de2 Mon Sep 17 00:00:00 2001 From: jefalon Date: Tue, 9 Jul 2024 13:52:33 -0600 Subject: [PATCH 06/15] fixed the double count of turns for reconstructing bcs --- tests/9-Regression/3D_Skew_Steady.yaml | 2 +- windse/DomainManager.py | 15 ++------------- 2 files changed, 3 insertions(+), 14 deletions(-) diff --git a/tests/9-Regression/3D_Skew_Steady.yaml b/tests/9-Regression/3D_Skew_Steady.yaml index 5f539d7d..472f3669 100644 --- a/tests/9-Regression/3D_Skew_Steady.yaml +++ b/tests/9-Regression/3D_Skew_Steady.yaml @@ -64,7 +64,7 @@ boundary_conditions: vel_profile: log HH_vel: 8.0 k: 0.4 - inflow_angle: 205.256 + inflow_angle: 205.25576915021696 problem: type: stabilized diff --git a/windse/DomainManager.py b/windse/DomainManager.py index b32c2d9b..416fe435 100644 --- a/windse/DomainManager.py +++ b/windse/DomainManager.py @@ -1126,23 +1126,12 @@ def RecomputeBoundaryMarkers(self,inflow_angle): ### Count the number of pi/2 sections in the new inflow_angle ### turns = inflow_angle/(pi/2) - # ### New order ### - # tol = 1e-3 - # if turns % 1 <= tol: # we are at a cardinal direction - # new_order = np.roll(cardinal_ids,int(turns)) - # else: - # new_order = np.roll(diagonal_ids,int(turns)) - ### New order ### tol = 1e-3 if turns % 1 <= tol: # we are at a cardinal direction - new_order = np.roll(cardinal_ids,int(np.sign(turns)*np.ceil(abs(turns)))) + new_order = np.roll(cardinal_ids,int(turns)) else: - new_order = np.roll(diagonal_ids,int(np.sign(turns)*np.ceil(abs(turns)))) - - - - + new_order = np.roll(diagonal_ids,int(turns)) ### Get a list of all wall facets ### wall_facets = [] From 7ae709fa13fefa422b9efc9e92dcf98da62725a2 Mon Sep 17 00:00:00 2001 From: jefalon Date: Tue, 9 Jul 2024 15:51:01 -0600 Subject: [PATCH 07/15] vectorized the angle clamp --- windse/helper_functions.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/windse/helper_functions.py b/windse/helper_functions.py index cf3253fa..7a077dea 100644 --- a/windse/helper_functions.py +++ b/windse/helper_functions.py @@ -21,10 +21,10 @@ def meteor_to_math(angle): ''' new_angle =np.radians(270.0-angle) - if new_angle<0: - new_angle+=2*np.pi - if new_angle>2*np.pi: - new_angle-=2*np.pi + # Clamp between 0 and 2pi + new_angle = (new_angle < 0) * (new_angle + 2*np.pi) \ + + np.logical_and(new_angle >= 0, new_angle <= 2*np.pi) * (new_angle) \ + + (new_angle > 2*np.pi) * (new_angle - 2*np.pi) return new_angle @@ -40,10 +40,10 @@ def math_to_meteor(angle): ''' new_angle =270.0-np.degrees(angle) - if new_angle<0: - new_angle+=360 - if new_angle>360: - new_angle-=360 + # Clamp between 0 and 360 + new_angle = (new_angle < 0) * (new_angle + 360) \ + + np.logical_and(new_angle >= 0, new_angle <= 360) * (new_angle) \ + + (new_angle > 360) * (new_angle - 360) return new_angle From 4ad85515c7dc174b68642145b537e9f3871c6e84 Mon Sep 17 00:00:00 2001 From: jefalon Date: Tue, 29 Oct 2024 21:43:13 -0600 Subject: [PATCH 08/15] added mock block --- windse/BoundaryManager.py | 99 ++++++++++++++++++++++- windse/ParameterManager.py | 7 +- windse/ProblemManager.py | 90 ++++++++++++++++++++- windse/RefinementManager.py | 10 ++- windse/SolverManager.py | 13 +-- windse/__init__.py | 4 +- windse/input_schema.yaml | 43 +++++++++- windse/wind_farm_types/GenericWindFarm.py | 7 +- windse/wind_farm_types/GridWindFarm.py | 3 +- windse/wind_farm_types/__init__.py | 1 + windse_driver/driver_functions.py | 1 + 11 files changed, 253 insertions(+), 25 deletions(-) diff --git a/windse/BoundaryManager.py b/windse/BoundaryManager.py index 7256fa49..9adbb8cd 100644 --- a/windse/BoundaryManager.py +++ b/windse/BoundaryManager.py @@ -79,6 +79,92 @@ def __init__(self,dom,fs,farm): for key,value in self.boundary_types.items(): self.boundary_types[key] = [i for i in value if i not in self.dom.bcs_to_remove] + self.u_bk = Function(self.fs.V) + self.bf_bk = Function(self.fs.V) + self.u_bg = Function(self.fs.V) + if self.use_bk_force == True: + print("doing this thing") + x = SpatialCoordinate(self.dom.mesh) + + # get the location of the neighboring farm + farm_RD = 130 + farm_offset_x = farm_RD*24 + farm_offset_y = farm_RD*6 + + # setup coeffs for either a body force or velocity + # bk_func_type = "velocity_perturbation" + # bk_func_type = "body_force" + if self.bk_func_type == "body_force": + gauss_coeffs = [ + [ + [0, 0, 6, 6, 0.03, 1.0, 1.0, 1.0], + ], + ] + elif self.bk_func_type == "velocity_perturbation": + gauss_coeffs = [ + [ + [-11, 0, 17, 22, -0.87, 1.0, 1.0, 0.25], + [ 8, -10, 25, 30, 0.49, 1.0, 1.0, 0.3], + [ 8, 10, 25, 30, 0.49, 1.0, 1.0, 0.3], + [ 30, 0, 35, 7, -1.0, 4.0, 4.0, 1.0], + ], + [ + [-3, 7, 17, 15, 0.6, 1.0, 1.0, 0.5], + [-3, -7, 17, 15, -0.6, 1.0, 1.0, 0.5], + [26, 7, 25, 11, -0.3, 1.0, 1.0, 0.25], + [26, -7, 25, 11, 0.3, 1.0, 1.0, 0.25], + ], + ] + # gauss_coeffs = [ + # [ + # [ -6, 0, 7, 9, -0.1, 1.0, 1.0, 1.0], + # [ 9, -20, 12, 12, 0.1, 1.0, 1.0, 1.0], + # [ 9, 20, 12, 12, 0.1, 1.0, 1.0, 1.0], + # [ 15, 0, 15, 9, -0.1, 3.0, 3.0, 1.0], + # ], + # [ + # [10, 7, 15, 15, 1.0, 1.0, 1.0, 1.0], + # [10, -7, 15, 15, -1.0, 1.0, 1.0, 1.0], + # [18, 7, 15, 15, -1.0, 1.0, 1.0, 1.0], + # [18, -7, 15, 15, 1.0, 1.0, 1.0, 1.0], + # ], + # ] + else: + raise ValueError(f"unknown type of blockage approximation: {self.bk_func_type}") + + # loop over coeff to build gauss functions + bk_func = [0,0,0] + for i in range(len(gauss_coeffs)): + comp_coeffs = gauss_coeffs[i] + + for j in range(len(comp_coeffs)): + val = comp_coeffs[j] + x0 = farm_RD*val[0] + farm_offset_x + y0 = farm_RD*val[1] + farm_offset_y + xs = farm_RD*val[2] + ys = farm_RD*val[3] + am = val[4] + px = val[5] + py = val[6] + po = val[7] + bk_func[i] += 1.0*am*exp(-( (((x[0]-x0)/xs)**2)**px + (((x[1]-y0)/ys)**2)**py )**po ) + + bk_func[i] *= exp(-((x[2]-90)/63)**2) + + # project the force to the velocity function space + bk_func = as_vector(bk_func) + + if self.bk_func_type == "body_force": + self.bf_bk = project(-bk_func,self.fs.V,solver_type='gmres',preconditioner_type="hypre_amg") + trash = self.params.Save(self.bf_bk,"bk_func",subfolder="functions/") + elif self.bk_func_type == "velocity_perturbation": + self.u_bk = project(-bk_func,self.fs.V,solver_type='gmres',preconditioner_type="hypre_amg") + trash = self.params.Save(self.u_bk,"bk_func",subfolder="functions/") + + + + + @no_annotations def DebugOutput(self): @@ -157,9 +243,8 @@ def SetupBoundaries(self): norm_comp = 1 elif bc_loc == 'bottom' or bc_loc == 'top': norm_comp = 2 - - self.bcu.append(DirichletBC(self.fs.V.sub(norm_comp), Constant(0.0), bc_domain)) - self.bcs.append(DirichletBC(self.fs.W.sub(0).sub(norm_comp), Constant(0.0), bc_domain)) + self.bcu.append(DirichletBC(self.fs.V.sub(norm_comp), self.u_bk, bc_domain)) + self.bcs.append(DirichletBC(self.fs.W.sub(0).sub(norm_comp), self.u_bk, bc_domain)) elif bc_type == 'no_stress': self.bcp.append(DirichletBC(self.fs.Q, Constant(0.0), bc_domain)) @@ -200,7 +285,9 @@ def SetupBoundaries(self): facet_normal = test_facet.normal().array() field_id = int(np.argmin(abs(abs(facet_normal)-1.0))) - bcu_eqns.append([self.fs.V.sub(field_id), self.fs.W.sub(0).sub(field_id), self.zero, boundary_id]) + temp = self.u_bk.split(deepcopy=True)[field_id] + temp.vector()[:] *= -0.5 + bcu_eqns.append([self.fs.V.sub(field_id), self.fs.W.sub(0).sub(field_id), temp, boundary_id]) elif bc_type == "no_stress": for b in bs: @@ -475,6 +562,10 @@ def __init__(self,dom,fs,farm): else: self.fs.VelocityAssigner.assign(self.bc_velocity,[self.ux,self.uy]) + ### Adjust by blockage function ### + self.bc_velocity.vector()[:] -= self.u_bk.vector()[:] + + ### Create Pressure Boundary Function self.bc_pressure = Function(self.fs.Q) diff --git a/windse/ParameterManager.py b/windse/ParameterManager.py index a0ef553c..c323d375 100644 --- a/windse/ParameterManager.py +++ b/windse/ParameterManager.py @@ -136,7 +136,10 @@ def TerminalUpdate(self,dic,keys,value): self.TerminalUpdate(next_dic,keys[1:],value) elif len(keys) == 1: current_value = dic.get(keys[0],"") - if isinstance(current_value,int): + # print(keys,value,current_value) + if isinstance(current_value,bool): + dic[keys[0]] = value == 'True' + elif isinstance(current_value,int): dic[keys[0]] = int(value) elif isinstance(current_value,float): dic[keys[0]] = float(value) @@ -148,7 +151,7 @@ def TerminalUpdate(self,dic,keys,value): if len(formatted_value) == 1 and "," in formatted_value[0]: formatted_value = formatted_value[0].split(",") - dic[keys[0]] = value + dic[keys[0]] = formatted_value def CheckParameters(self,updates,defaults,out_string=""): diff --git a/windse/ProblemManager.py b/windse/ProblemManager.py index 4fa6f398..7abe2f86 100644 --- a/windse/ProblemManager.py +++ b/windse/ProblemManager.py @@ -352,12 +352,13 @@ def ComputeFunctional(self,inflow_angle): self.ReyStress=self.nu_T*grad(self.u_k) self.vertKE= self.ReyStress[0,2]*self.u_k[0] + ### Create the functional ### - self.F = inner(grad(self.u_k)*self.u_k, v)*dx - self.F += Sx*Sx*(nu+self.nu_T)*inner(grad(self.u_k), grad(v))*dx + self.F = inner(grad(self.u_k+self.bd.u_bk)*(self.u_k+self.bd.u_bk), v)*dx + self.F += Sx*Sx*(nu+self.nu_T)*inner(grad(self.u_k+self.bd.u_bk), grad(v))*dx self.F += - inner(div(v),self.p_k)*dx - self.F += - inner(div(self.u_k),q)*dx - self.F += - inner(f,v)*dx + self.F += - inner(div(self.u_k+self.bd.u_bk),q)*dx + self.F += - inner(self.bd.bf_bk,v)*dx self.F += - tf_term @@ -366,6 +367,85 @@ def ComputeFunctional(self,inflow_angle): self.fprint("Using Body Force") self.F += inner(-self.mbody_force*self.bd.inflow_unit_vector,v)*dx + # print(self.use_bk_force) + # if self.use_bk_force == True: + # print("doing this thing") + # x = SpatialCoordinate(self.dom.mesh) + + # # x_test_range = np.linspace(-378, 2268,10000) + # # y_test_range = np.linspace(-1890, 1890,10000) + + # # # bg_force_S = (x[0]/1500*x[1]/2000)**2.0 + # # # self.F += -inner(div(v),bg_force_S)*dx + # # # self.F += inner(v,grad(bg_force_S))*dx + # # # self.F += inner(v,bg_force_V)*dx + + # # bg_force_V = as_tensor(((x[0]/1500*x[1]/2000)*1/1500*x[1]/2000,(x[0]/1500*x[1]/2000)*x[0]/1500*1/2000,0.0)) + # # self.F += inner(v,bg_force_V)*dx + + # # x_0 = -1500.0 + # # x_1 = -1500.0 + # # x_mag = 1.0/max((x_test_range-x_0)*(x_test_range-x_1)*(x_test_range-x_0)*(x_test_range-x_1)) + # # y_0 = -1000.0 + # # y_1 = 1000.0 + # # y_mag = 1.0/max((y_test_range-y_0)*(y_test_range-y_1)) + # # bg_force = -100.0* x_mag*(x[0]-x_0)*(x[0]-x_1)*(x[0]-x_0)*(x[0]-x_1) * y_mag*(x[1]-y_0)*(x[1]-y_1) + + # # bg_force_save = project(bg_force,self.fs.Q,solver_type='gmres',preconditioner_type="hypre_amg") + # # trash = self.params.Save(bg_force_save,"bg_force",subfolder="functions/") + + # # # self.F += inner(bg_force,q)*dx + # # self.F += -inner(div(v),bg_force)*dx + + + # x_test_range = np.linspace(-3780, 2268,10000) + # y_test_range = np.linspace(-3780, 3780,10000) + + # x_0 = -3780. + # x_1 = -3780. + # x_mag = 1.0/max((x_test_range-x_0)*(x_test_range-x_1)) + # y_0 = -1764.0+756 + # y_1 = 1764.0+756 + # y_mag = 1.0/max((y_test_range-y_0)*(y_test_range-y_1))#*(y_test_range-y_2)*(y_test_range-y_3))) + # x_force = 0.0005* (x_mag*(x[0]-x_0)*(x[0]-x_1)) * (y_mag*(x[1]-y_0)*(x[1]-y_1))#*(x[1]-y_2)*(x[1]-y_3)) + + # x_0 = -3780. + # x_1 = -3780. + # x_mag = 1.0/max((x_test_range-x_0)*(x_test_range-x_1)) + # y_0 = -3780.0+756 + # y_1 = 0.0+756 + # y_2 = 3780.0+756 + # y_mag = 1.0/max((y_test_range-y_0)*(y_test_range-y_1)*(y_test_range-y_2)) + # y_force = 0.00005* (x_mag*(x[0]-x_0)*(x[0]-x_1))**2.0 * y_mag*(x[1]-y_0)*(x[1]-y_1)*(x[1]-y_2) + + # # # x_W = 700.0 + # # # y_W = 700.0 + # # # x_0 = 0.0 + # # # y_0 = 0.0 + # # # x_force = 0.05*exp(-pow(((x[0]-x_0)/x_W)**2+((x[1]-y_0)/y_W)**2,1.0)) + + + # z_force = 0.0 + # bg_force = as_vector((x_force,y_force,z_force)) + + # bg_force_save = project(bg_force,self.fs.V,solver_type='gmres',preconditioner_type="hypre_amg") + # trash = self.params.Save(bg_force_save,"bg_force",subfolder="functions/") + + + # # # nu_T_background=self.ComputeTurbulenceModel(bg_force) + # # # F_background = 0 + # # # F_background += inner(grad(bg_force)*bg_force, v)*dx + # # # F_background += Sx*Sx*(nu+nu_T_background)*inner(grad(bg_force), grad(v))*dx + + + + + + + # # # self.F += -F_background + # self.F += inner(-bg_force,v)*dx + + ################ THIS IS A CHEAT ####################\ if self.use_corrective_force: self.fprint("Using Corrective Force") @@ -403,6 +483,8 @@ def ComputeFunctional(self,inflow_angle): else: term25 = (abs(sin(inflow_angle))*dudx*q + abs(cos(inflow_angle))*dvdy*q)*dx + + self.F -= term25 self.fprint("Stabilized Problem Setup",special="footer") diff --git a/windse/RefinementManager.py b/windse/RefinementManager.py index 5384c83f..ce9433bd 100644 --- a/windse/RefinementManager.py +++ b/windse/RefinementManager.py @@ -35,8 +35,8 @@ def CreateRefinementList(dom, farm, refine_params): for key in sorted_refine_custom_keys: refine_list.append(refine_custom[key]) - if not keys_are_integers: - refine_list[-1]["type"] = key + # if not keys_are_integers: + # refine_list[-1]["type"] = key if farm_num > 0: bbox = farm.calculate_farm_bounding_box() @@ -95,7 +95,8 @@ def CreateRefinementList(dom, farm, refine_params): "theta": theta, "pivot_offset": pivot_offset, "expand_factor": expand_factor}) - + else: + raise ValueError(f"unknown farm refine: {farm_type}") if turbine_num > 0 and farm.numturbs > 0: for i in range(turbine_num,0,-1): expand_factor = (turbine_factor)**(i) @@ -133,6 +134,8 @@ def CreateRefinementList(dom, farm, refine_params): "radius": radius, "theta": theta, "expand_factor": expand_factor}) + else: + raise ValueError(f"unknown turbine refine: {turbine_type}") if refine_power_calc: radius = max(RDs) @@ -147,7 +150,6 @@ def CreateRefinementList(dom, farm, refine_params): "expand_factor": expand_factor, "centered": centered}) - return refine_list def RefineMesh(dom,farm): diff --git a/windse/SolverManager.py b/windse/SolverManager.py index c76dff93..2c5c3e8a 100644 --- a/windse/SolverManager.py +++ b/windse/SolverManager.py @@ -78,6 +78,7 @@ def __init__(self,problem): self.simTime_prev = None self.iter_val = 0.0 self.pow_saved = False + self.save_extra = False ### Update attributes based on params file ### for key, value in self.params["solver"].items(): @@ -162,16 +163,16 @@ def Save(self,val=0): if self.first_save: self.u_file = self.params.Save(u,"velocity",subfolder="solutions/",val=val) self.p_file = self.params.Save(p,"pressure",subfolder="solutions/",val=val) - self.nuT_file = self.params.Save(self.nu_T,"eddy_viscosity",subfolder="solutions/",val=val) - if self.problem.dom.dim == 3: + if self.problem.dom.dim == 3 and self.save_extra: + self.nuT_file = self.params.Save(self.nu_T,"eddy_viscosity",subfolder="solutions/",val=val) self.ReyStress_file = self.params.Save(self.ReyStress,"Reynolds_stresses",subfolder="solutions/",val=val) self.vertKE_file = self.params.Save(self.vertKE,"Vertical KE",subfolder="solutions/",val=val) self.first_save = False else: self.params.Save(u,"velocity",subfolder="solutions/",val=val,file=self.u_file) self.params.Save(p,"pressure",subfolder="solutions/",val=val,file=self.p_file) - self.params.Save(self.nu_T,"eddy_viscosity",subfolder="solutions/",val=val,file=self.nuT_file) - if self.problem.dom.dim == 3: + if self.problem.dom.dim == 3 and self.save_extra: + self.params.Save(self.nu_T,"eddy_viscosity",subfolder="solutions/",val=val,file=self.nuT_file) self.params.Save(self.ReyStress,"Reynolds_stresses",subfolder="solutions/",val=val,file=self.ReyStress_file) self.params.Save(self.vertKE,"Vertical KE",subfolder="solutions/",val=val,file=self.vertKE_file) u.vector()[:]=u.vector()[:]*self.problem.dom.xscale @@ -458,10 +459,10 @@ def Solve(self): ### Hack into doflin adjoint to update the local controls at the start of the adjoint solve ### - self.nu_T = project(self.problem.nu_T,self.problem.fs.Q,solver_type='gmres',preconditioner_type="hypre_amg",**self.extra_kwarg) - if self.problem.dom.dim == 3: + if self.problem.dom.dim == 3 and self.save_extra: self.fprint("") self.fprint("Projecting Reynolds Stress") + self.nu_T = project(self.problem.nu_T,self.problem.fs.Q,solver_type='gmres',preconditioner_type="hypre_amg",**self.extra_kwarg) self.ReyStress = project(self.problem.ReyStress,self.problem.fs.T,solver_type='gmres',preconditioner_type="hypre_amg",**self.extra_kwarg) self.vertKE = project(self.problem.vertKE,self.problem.fs.Q,solver_type='gmres',preconditioner_type="hypre_amg",**self.extra_kwarg) diff --git a/windse/__init__.py b/windse/__init__.py index 785bc911..72435a55 100644 --- a/windse/__init__.py +++ b/windse/__init__.py @@ -37,8 +37,8 @@ def initialize(loc,updated_parameters=[]): global BoxDomain, CylinderDomain, CircleDomain, RectangleDomain, ImportedDomain, InterpolatedCylinderDomain, InterpolatedBoxDomain, PeriodicDomain from windse.DomainManager import BoxDomain, CylinderDomain, CircleDomain, RectangleDomain, ImportedDomain, InterpolatedCylinderDomain, InterpolatedBoxDomain, PeriodicDomain - global GridWindFarm, RandomWindFarm, ImportedWindFarm, EmptyWindFarm - from windse.wind_farm_types import GridWindFarm, RandomWindFarm, ImportedWindFarm, EmptyWindFarm + global GridWindFarm, TrapWindFarm, RandomWindFarm, ImportedWindFarm, EmptyWindFarm + from windse.wind_farm_types import GridWindFarm, TrapWindFarm, RandomWindFarm, ImportedWindFarm, EmptyWindFarm global RefineMesh, WarpMesh from windse.RefinementManager import RefineMesh, WarpMesh diff --git a/windse/input_schema.yaml b/windse/input_schema.yaml index bb0356d6..7b94601b 100644 --- a/windse/input_schema.yaml +++ b/windse/input_schema.yaml @@ -344,6 +344,7 @@ properties: - "null" enum: - "grid" + - "trap" - "random" - "imported" - "empty" @@ -420,6 +421,32 @@ properties: - "null" description: "Integer number of turbines in the x direction." units: "dimensionless" + num_x: + default: Null + type: + - "integer" + - "null" + description: "Integer number of turbines in the x direction." + units: "dimensionless" + num_y: + default: Null + type: + - "integer" + - "null" + description: "Integer number of turbines in the y direction." + units: "dimensionless" + expand_factor: + default: 1.0 + type: + - "number" + description: "How much the farm grows (>1) or shrinks (<1) between the first and last columns." + units: "dimensionless" + leading_center: + default: [0.0,0.0] + type: + - "array" + description: "Location of the center of the leading column of turbines." + units: "meter" jitter: default: 0.0 type: @@ -648,7 +675,7 @@ properties: description: "Integer number of farm level refinements." units: "dimensionless" farm_type: - default: "square" + default: "box" type: "string" enum: - "full" @@ -880,6 +907,20 @@ properties: default: False type: "boolean" description: "Relax divergence free constraint to entrain momentum from 'above' and 'below'." + use_bk_force: + default: False + type: "boolean" + description: "When True, adds a background body force to mimic blockage of near by wind farm." + bk_func_type: + default: Null + type: + - "string" + - "null" + description: "choose the type of blockage approximation." + enum: + - "body_force" + - "velocity_perturbation" + - Null viscosity: default: 0.1 type: "number" diff --git a/windse/wind_farm_types/GenericWindFarm.py b/windse/wind_farm_types/GenericWindFarm.py index 3b412879..124cc4c3 100644 --- a/windse/wind_farm_types/GenericWindFarm.py +++ b/windse/wind_farm_types/GenericWindFarm.py @@ -2,7 +2,7 @@ from windse.turbine_types import turbine_dict import numpy as np import time, os -from . import MeshFunction, CompiledSubDomain, Measure, cells, project, inner, FiniteElement, FunctionSpace, MixedElement, assemble, dx, parameters, Form, File +from . import MeshFunction, CompiledSubDomain, Measure, cells, project, inner, FiniteElement, FunctionSpace, MixedElement, assemble, dx, parameters, Form, File, Function import matplotlib.pyplot as plt from pyadjoint.tape import stop_annotating from pyadjoint import AdjFloat @@ -199,6 +199,10 @@ def compute_turbine_force(self,u,v,inflow_angle,fs,**kwargs): # add the subdomain to the marker subdomain_marker.mark(self.turbine_subdomains,turb.index+1) + # handle the case where there are no turbines + if self.numturbs == 0: + print("testing") + self.tf_list = [Function(self.fs.V),] # create a measure for the turbines self.local_dx = Measure('dx', subdomain_data=self.turbine_subdomains) @@ -218,6 +222,7 @@ def compute_turbine_force(self,u,v,inflow_angle,fs,**kwargs): # use the full domain of integration if local is not requested else: + # print(sum(final_tf_list)) tf_term = inner(sum(final_tf_list),v)*dx return tf_term diff --git a/windse/wind_farm_types/GridWindFarm.py b/windse/wind_farm_types/GridWindFarm.py index 4a954f6a..8b1b07cb 100644 --- a/windse/wind_farm_types/GridWindFarm.py +++ b/windse/wind_farm_types/GridWindFarm.py @@ -59,7 +59,8 @@ def compute_parameters(self): if self.ex_x is None: x_dist = self.x_spacing * (self.grid_cols - 1) y_dist = self.y_spacing * (self.grid_rows - 1) - self.ex_x = [0., x_dist + 2 * self.radius] + # self.ex_x = [0., x_dist + 2 * self.radius] + self.ex_x = [-x_dist / 2 - self.radius, x_dist / 2 + self.radius] self.ex_y = [-y_dist / 2 - self.radius, y_dist / 2 + self.radius] else: self.ex_x = [self.ex_x[0] + self.radius, self.ex_x[1] - self.radius] diff --git a/windse/wind_farm_types/__init__.py b/windse/wind_farm_types/__init__.py index a508f8c1..c9b5f98a 100644 --- a/windse/wind_farm_types/__init__.py +++ b/windse/wind_farm_types/__init__.py @@ -10,6 +10,7 @@ ### import the wind farm types from .GenericWindFarm import GenericWindFarm from .GridWindFarm import GridWindFarm +from .TrapWindFarm import TrapWindFarm from .RandomWindFarm import RandomWindFarm from .ImportedWindFarm import ImportedWindFarm from .EmptyWindFarm import EmptyWindFarm diff --git a/windse_driver/driver_functions.py b/windse_driver/driver_functions.py index cd66da3c..acedb2e5 100644 --- a/windse_driver/driver_functions.py +++ b/windse_driver/driver_functions.py @@ -109,6 +109,7 @@ def BuildDomain(params): #### Build Farm farm_dict = {"grid":windse.GridWindFarm, + "trap":windse.TrapWindFarm, "random":windse.RandomWindFarm, "imported":windse.ImportedWindFarm, "empty":windse.EmptyWindFarm} From 1d590d250dbcd02f5e146b62db19b3104f8b8264 Mon Sep 17 00:00:00 2001 From: jefalon Date: Tue, 29 Oct 2024 22:38:53 -0600 Subject: [PATCH 09/15] added trap farm --- windse/wind_farm_types/TrapWindFarm.py | 87 ++++++++++++++++++++++++++ 1 file changed, 87 insertions(+) create mode 100644 windse/wind_farm_types/TrapWindFarm.py diff --git a/windse/wind_farm_types/TrapWindFarm.py b/windse/wind_farm_types/TrapWindFarm.py new file mode 100644 index 00000000..e170aba8 --- /dev/null +++ b/windse/wind_farm_types/TrapWindFarm.py @@ -0,0 +1,87 @@ +from . import GenericWindFarm +import numpy as np + +class TrapWindFarm(GenericWindFarm): + """ + A TrapWindFarm produces turbines on a grid but shape as a trapezoid. The params.yaml file determines + how this grid is set up. + + Example: + In the .yaml file you need to define:: + + wind_farm: + # # Description | Units + type: trap + num_x: # Number of columns | - + num_y: # Number of rows | - + x_spacing: 300 # Space between columns | m + y_spacing: 300 # Space between row | m + expand_factor: 2 # How much the farm grows (>1) or shrinks (<1) between the first and last columns + leading_center:[0,0] # location of the center of the leading column of turbines + + + This will produce a 6x6 grid of turbines where the first column is + 300 meters between turbines, and the last column is 600 m between + turbines. The first column will be centered on the point 0,0 + + Args: + dom (:meth:`windse.DomainManager.GenericDomain`): a windse domain object. + """ + def __init__(self,dom): + + self.name = "Trap Farm" + super(TrapWindFarm, self).__init__(dom) + + def load_parameters(self): + + # Load from the yaml file + self.num_x = self.params["wind_farm"]["num_x"] + self.num_y = self.params["wind_farm"]["num_y"] + self.x_spacing = self.params["wind_farm"]["x_spacing"] + self.y_spacing = self.params["wind_farm"]["y_spacing"] + self.expand_factor = self.params["wind_farm"]["expand_factor"] + self.leading_center = self.params["wind_farm"]["leading_center"] + + def compute_parameters(self): + + # compute number of turbines + self.numturbs = self.num_x*self.num_y + + def initialize_turbine_locations(self): + + ### Create the x and y coords ### + if self.num_x > 1: + x_list = np.linspace(0,1,self.num_x) + else: + raise ValueError("the number of column for a TrapWindFarm needs to be greater than 1") + if self.num_y > 1: + y_list = np.linspace(-1,1,self.num_y) + else: + raise ValueError("the number of rows for a TrapWindFarm needs to be greater than 1") + + ### Use the x and y coords to make a mesh grid ### + X, Y = np.meshgrid(x_list,y_list) + + # Distort grid such that the first col remains unchanged but the last one has the right spacing + Y = Y*(self.expand_factor*X-X+1) + + # Scale col to real locations + X = self.x_spacing*(self.num_x-1)*(X)+self.leading_center[0] + + # Scale row to real locations + Y = self.y_spacing*(self.num_y-1)*(Y/2)+self.leading_center[1] + + # Convert to lists + x = X.flatten() + y = Y.flatten() + + # Recompute the extents based on grid manipulations + self.ex_x = [np.min(x), np.max(x)] + self.ex_y = [np.min(y), np.max(y)] + + # Output some useful data + self.fprint("Grid Size: {:d} x {:d}".format(self.num_y,self.num_x)) + self.fprint("X Range: [{: 1.2f}, {: 1.2f}]".format(self.ex_x[0],self.ex_x[1])) + self.fprint("Y Range: [{: 1.2f}, {: 1.2f}]".format(self.ex_y[0],self.ex_y[1])) + + return np.array([x,y]).T \ No newline at end of file From a0aacda2608a7fe07d04b8a450a756418aedcb39 Mon Sep 17 00:00:00 2001 From: jefalon Date: Sun, 3 Nov 2024 16:37:10 -0700 Subject: [PATCH 10/15] made more neighbor farm yaml values --- windse/BoundaryManager.py | 29 ++++++++++++----------- windse/SolverManager.py | 1 + windse/input_schema.yaml | 9 +++++++ windse/wind_farm_types/GenericWindFarm.py | 1 + 4 files changed, 26 insertions(+), 14 deletions(-) diff --git a/windse/BoundaryManager.py b/windse/BoundaryManager.py index 9adbb8cd..e2bbbae8 100644 --- a/windse/BoundaryManager.py +++ b/windse/BoundaryManager.py @@ -83,13 +83,14 @@ def __init__(self,dom,fs,farm): self.bf_bk = Function(self.fs.V) self.u_bg = Function(self.fs.V) if self.use_bk_force == True: - print("doing this thing") + self.fprint(f"Applying Mock Block: {self.bk_func_type}") + self.fprint(f"Neighbor Offset: {self.neighbor_offset}") x = SpatialCoordinate(self.dom.mesh) # get the location of the neighboring farm - farm_RD = 130 - farm_offset_x = farm_RD*24 - farm_offset_y = farm_RD*6 + farm_RD = np.mean(farm.get_rotor_diameters()) + farm_offset_x = farm_RD*self.neighbor_offset[0] + farm_offset_y = farm_RD*self.neighbor_offset[1] # setup coeffs for either a body force or velocity # bk_func_type = "velocity_perturbation" @@ -97,22 +98,22 @@ def __init__(self,dom,fs,farm): if self.bk_func_type == "body_force": gauss_coeffs = [ [ - [0, 0, 6, 6, 0.03, 1.0, 1.0, 1.0], + [0, 0, 6, 6, 0.04, 1.0, 1.0, 1.0], ], ] elif self.bk_func_type == "velocity_perturbation": gauss_coeffs = [ [ - [-11, 0, 17, 22, -0.87, 1.0, 1.0, 0.25], - [ 8, -10, 25, 30, 0.49, 1.0, 1.0, 0.3], - [ 8, 10, 25, 30, 0.49, 1.0, 1.0, 0.3], - [ 30, 0, 35, 7, -1.0, 4.0, 4.0, 1.0], + [-11, 0, 17, 22, -0.87*1.5, 1.0, 1.0, 0.25], + [ 8, -10, 25, 30, 0.49*1.5, 1.0, 1.0, 0.3], + [ 8, 10, 25, 30, 0.49*1.5, 1.0, 1.0, 0.3], + [ 30, 0, 35, 7, -1.0*1.5, 4.0, 4.0, 1.0], ], [ - [-3, 7, 17, 15, 0.6, 1.0, 1.0, 0.5], - [-3, -7, 17, 15, -0.6, 1.0, 1.0, 0.5], - [26, 7, 25, 11, -0.3, 1.0, 1.0, 0.25], - [26, -7, 25, 11, 0.3, 1.0, 1.0, 0.25], + [-3, 7, 17, 15, 0.6*1.5, 1.0, 1.0, 0.5], + [-3, -7, 17, 15, -0.6*1.5, 1.0, 1.0, 0.5], + [26, 7, 25, 11, -0.3*1.5, 1.0, 1.0, 0.25], + [26, -7, 25, 11, 0.3*1.5, 1.0, 1.0, 0.25], ], ] # gauss_coeffs = [ @@ -149,7 +150,7 @@ def __init__(self,dom,fs,farm): po = val[7] bk_func[i] += 1.0*am*exp(-( (((x[0]-x0)/xs)**2)**px + (((x[1]-y0)/ys)**2)**py )**po ) - bk_func[i] *= exp(-((x[2]-90)/63)**2) + bk_func[i] *= exp(-((x[2]-self.vel_height)/(farm_RD/2))**2) # project the force to the velocity function space bk_func = as_vector(bk_func) diff --git a/windse/SolverManager.py b/windse/SolverManager.py index 2c5c3e8a..070f346d 100644 --- a/windse/SolverManager.py +++ b/windse/SolverManager.py @@ -229,6 +229,7 @@ def EvaulatePowerFunctional(self): "simTime": self.simTime } out = self.problem.farm.save_power(self.problem.u_k,self.problem.dom.inflow_angle, **kwargs) + self.fprint(f"Total Power: {out}") return out diff --git a/windse/input_schema.yaml b/windse/input_schema.yaml index 7b94601b..22e8bf73 100644 --- a/windse/input_schema.yaml +++ b/windse/input_schema.yaml @@ -921,6 +921,15 @@ properties: - "body_force" - "velocity_perturbation" - Null + neighbor_offset: + default: [0,0] + type: + - "object" + description: "location in rotor diameters of the neighboring farm" + enum: + - "body_force" + - "velocity_perturbation" + - Null viscosity: default: 0.1 type: "number" diff --git a/windse/wind_farm_types/GenericWindFarm.py b/windse/wind_farm_types/GenericWindFarm.py index 124cc4c3..4041c3e3 100644 --- a/windse/wind_farm_types/GenericWindFarm.py +++ b/windse/wind_farm_types/GenericWindFarm.py @@ -425,6 +425,7 @@ def save_power(self, u, inflow_angle, iter_val = 0.0, simTime = 0.0): self.power_first_save = False else: self.params.save_csv("power_data",data=[J_list],subfolder=self.params.folder+"data/",mode='a') + return J_list[-1] def debug_output(self): """ From 0e600d32ae47d114f07236ca812f8acb5114a7a1 Mon Sep 17 00:00:00 2001 From: jefalon Date: Tue, 1 Apr 2025 11:08:19 -0600 Subject: [PATCH 11/15] update blockage force definitions --- windse/BoundaryManager.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/windse/BoundaryManager.py b/windse/BoundaryManager.py index e2bbbae8..d1b14856 100644 --- a/windse/BoundaryManager.py +++ b/windse/BoundaryManager.py @@ -98,22 +98,22 @@ def __init__(self,dom,fs,farm): if self.bk_func_type == "body_force": gauss_coeffs = [ [ - [0, 0, 6, 6, 0.04, 1.0, 1.0, 1.0], + [0, 0, 6, 6, 0.02, 1.0, 1.0, 1.0], ], ] elif self.bk_func_type == "velocity_perturbation": gauss_coeffs = [ [ - [-11, 0, 17, 22, -0.87*1.5, 1.0, 1.0, 0.25], - [ 8, -10, 25, 30, 0.49*1.5, 1.0, 1.0, 0.3], - [ 8, 10, 25, 30, 0.49*1.5, 1.0, 1.0, 0.3], - [ 30, 0, 35, 7, -1.0*1.5, 4.0, 4.0, 1.0], + [-11, 0, 17, 22, -0.87, 1.0, 1.0, 0.25], + [ 8, -10, 25, 30, 0.49, 1.0, 1.0, 0.3], + [ 8, 10, 25, 30, 0.49, 1.0, 1.0, 0.3], + [ 30, 0, 35, 7, -1.0, 4.0, 4.0, 1.0], ], [ - [-3, 7, 17, 15, 0.6*1.5, 1.0, 1.0, 0.5], - [-3, -7, 17, 15, -0.6*1.5, 1.0, 1.0, 0.5], - [26, 7, 25, 11, -0.3*1.5, 1.0, 1.0, 0.25], - [26, -7, 25, 11, 0.3*1.5, 1.0, 1.0, 0.25], + [-3, 7, 17, 15, 0.6, 1.0, 1.0, 0.5], + [-3, -7, 17, 15, -0.6, 1.0, 1.0, 0.5], + [26, 7, 25, 11, -0.3, 1.0, 1.0, 0.25], + [26, -7, 25, 11, 0.3, 1.0, 1.0, 0.25], ], ] # gauss_coeffs = [ From c2d5390d2624f5bcb0feaa799da9a8c5984d1bbe Mon Sep 17 00:00:00 2001 From: jefalon Date: Tue, 29 Apr 2025 12:37:41 -0600 Subject: [PATCH 12/15] added neighbor farm approximation yaml options --- windse/input_schema.yaml | 42 ++++++++++++++++++---------------------- 1 file changed, 19 insertions(+), 23 deletions(-) diff --git a/windse/input_schema.yaml b/windse/input_schema.yaml index 22e8bf73..63bca1ec 100644 --- a/windse/input_schema.yaml +++ b/windse/input_schema.yaml @@ -882,6 +882,25 @@ properties: - "array" - "null" description: "List of no_stress bc names." + use_bk_force: + default: False + type: "boolean" + description: "When True, adds a background body force to mimic blockage of near by wind farm." + bk_func_type: + default: Null + type: + - "string" + - "null" + description: "choose the type of blockage approximation." + enum: + - "body_force" + - "velocity_perturbation" + - Null + neighbor_offset: + default: [0,0] + type: + - "object" + description: "location in rotor diameters of the neighboring farm" # ================================================================ problem: # additionalProperties: False @@ -907,29 +926,6 @@ properties: default: False type: "boolean" description: "Relax divergence free constraint to entrain momentum from 'above' and 'below'." - use_bk_force: - default: False - type: "boolean" - description: "When True, adds a background body force to mimic blockage of near by wind farm." - bk_func_type: - default: Null - type: - - "string" - - "null" - description: "choose the type of blockage approximation." - enum: - - "body_force" - - "velocity_perturbation" - - Null - neighbor_offset: - default: [0,0] - type: - - "object" - description: "location in rotor diameters of the neighboring farm" - enum: - - "body_force" - - "velocity_perturbation" - - Null viscosity: default: 0.1 type: "number" From 0474252fb9ef10ec5f75311d8f3b3cf1241006e0 Mon Sep 17 00:00:00 2001 From: agushin101 Date: Tue, 6 May 2025 18:39:15 -0400 Subject: [PATCH 13/15] build new environment --- environment.yaml | 14 +++++++------- install.sh | 32 ++++++++++++++++++++++---------- 2 files changed, 29 insertions(+), 17 deletions(-) diff --git a/environment.yaml b/environment.yaml index 8fc0ec58..433c3548 100644 --- a/environment.yaml +++ b/environment.yaml @@ -3,17 +3,15 @@ channels: - conda-forge - defaults dependencies: - - python=3.9 #3.10 causes CI to hang in parallel - - mpich - # - mpich=4.1.2=hd33e60e_101 - # - mpich=4.0.3 + # - python=3.7 - coveralls - dolfin-adjoint=2019.1.0 - - fenics + - fenics=2019.1.0=py39*_26 - hdf5 - jsonschema - matplotlib - memory_profiler + - mshr - openmdao - pandas - pyoptsparse @@ -21,7 +19,8 @@ dependencies: - pytest-cov - pytest-mpi - pyyaml - - scipy + - scipy=1.7.0 + - sympy=1.13.3 - slepc - sphinx<7.0.0 - pip @@ -37,5 +36,6 @@ dependencies: - pulp - fatpack - gmsh - - meshio + - meshio==5.3.4 + - h5py==3.13.0 - -e . diff --git a/install.sh b/install.sh index 1b2457f2..a22cf661 100644 --- a/install.sh +++ b/install.sh @@ -1,15 +1,27 @@ #!/bin/bash -### Create some names -conda_base_path=$(conda info --base) -env_name=$1 - ### Run this to use conda in the script -source ${conda_base_path}/etc/profile.d/conda.sh +source $(conda info --base)/etc/profile.d/conda.sh + +# ### Create the Environment +# conda create -y --name $1 +# conda activate $1 + +# ### Install conda-forge dependencies +# conda install -y -c conda-forge fenics=2019.1.0=py38_9 dolfin-adjoint matplotlib scipy=1.4.1 slepc mshr hdf5 pyyaml memory_profiler pytest pytest-cov pytest-mpi coveralls pandas + +# ### Install the new tsfc compiler +# pip install git+https://github.com/blechta/tsfc.git@2018.1.0 +# pip install git+https://github.com/blechta/COFFEE.git@2018.1.0 +# pip install git+https://github.com/blechta/FInAT.git@2018.1.0 +# pip install git+https://github.com/mdolab/pyoptsparse@v1.0 +# pip install singledispatch networkx pulp openmdao fatpack + +# ### Install editible version of WindSE +# pip install -e . + +conda env create --name ${1} --file environment.yaml -### Create the Environment with mamba -conda create --yes --name ${env_name} --channel conda-forge python=3.9 mamba -conda activate ${env_name} +conda activate ${1} -### Install dependencies from environment -mamba env update --prefix ${conda_base_path}/envs/${env_name} --file environment.yaml --prune +pip install aeromesh --no-deps \ No newline at end of file From 91131732eb04b02471e3c1a2a318d4f010b90fa6 Mon Sep 17 00:00:00 2001 From: jefalon Date: Fri, 9 May 2025 12:52:44 -0600 Subject: [PATCH 14/15] added a check for the 'windse run' command in parallel --- windse_driver/driver.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/windse_driver/driver.py b/windse_driver/driver.py index c492bd6c..85117361 100644 --- a/windse_driver/driver.py +++ b/windse_driver/driver.py @@ -1,5 +1,6 @@ import os import __main__ +from mpi4py import MPI ### Get the name of program importing this package ### if hasattr(__main__,"__file__"): @@ -54,6 +55,12 @@ def get_action(): ### Run the driver ### def run_action(params_loc=None): + + + # Check to see if call from exec and parallel + if __name__ == "windse_driver.driver" and MPI.COMM_WORLD.Get_size()>1: + raise NotImplementedError("Using the 'windse run' command in parallel is not currently supported. Instead run 'mpirun -n # python -u WindSE/windse_driver/driver.py run '") + tick = time.time() ### Clean up other module references ### From fd879387d8e32116c774805d8708b21c7bbae089 Mon Sep 17 00:00:00 2001 From: jefalon Date: Tue, 13 May 2025 13:36:25 -0600 Subject: [PATCH 15/15] removed mpi4py import --- windse_driver/driver.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/windse_driver/driver.py b/windse_driver/driver.py index 85117361..1b73a65f 100644 --- a/windse_driver/driver.py +++ b/windse_driver/driver.py @@ -1,6 +1,5 @@ import os import __main__ -from mpi4py import MPI ### Get the name of program importing this package ### if hasattr(__main__,"__file__"): @@ -58,7 +57,7 @@ def run_action(params_loc=None): # Check to see if call from exec and parallel - if __name__ == "windse_driver.driver" and MPI.COMM_WORLD.Get_size()>1: + if __name__ == "windse_driver.driver" and dolfin.MPI.comm_world.Get_size()>1: raise NotImplementedError("Using the 'windse run' command in parallel is not currently supported. Instead run 'mpirun -n # python -u WindSE/windse_driver/driver.py run '") tick = time.time()