diff --git a/environment.yaml b/environment.yaml index 8fc0ec5..433c354 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 1b2457f..a22cf66 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 diff --git a/tests/9-Regression/Truth_Data/tolerances.yaml b/tests/9-Regression/Truth_Data/tolerances.yaml index 6c4f2ad..698ff82 100644 --- a/tests/9-Regression/Truth_Data/tolerances.yaml +++ b/tests/9-Regression/Truth_Data/tolerances.yaml @@ -6,4 +6,6 @@ SolverManager: "*max": [1e-3,"absolute"] "*min": [1e-3,"absolute"] DomainManager: - volume: [1e-9, "relative"] \ No newline at end of file + volume: [1e-9, "relative"] +BoundaryManager: + avg_x: [1e-3, "relative"] \ No newline at end of file diff --git a/windse/BoundaryManager.py b/windse/BoundaryManager.py index 7256fa4..671b3b9 100644 --- a/windse/BoundaryManager.py +++ b/windse/BoundaryManager.py @@ -79,6 +79,104 @@ 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: + 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 = 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" + # bk_func_type = "body_force" + if self.bk_func_type == "body_force": + # gauss_coeffs = [ + # [ + # [0, 0, 6, 6, 0.02, 1.0, 1.0, 1.0], + # ], + gauss_coeffs = [ + [ + [0, -30, 40, 20, 0.0045, 1.0, 1.0, 1.0], + # [0, -31.5, 45, 23.25, 0.005, 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 + if self.dom.dim == 2: + bk_func = [0,0] + if self.dom.dim == 3: + 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 ) + + if self.dom.dim == 3: + 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) + + 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 +255,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 +297,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 +574,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 a0ef553..c323d37 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 4fa6f39..b9dc25b 100644 --- a/windse/ProblemManager.py +++ b/windse/ProblemManager.py @@ -232,6 +232,15 @@ def SimpleControlUpdate(self): self.u_k, self.p_k = split(self.up_k) self.farm.SimpleControlUpdate() + def ChangeAxial(self,axial): + adj_start = time.time() + self.fprint("Adjusting Axial Induction",special="header") + self.fprint("New Value: {:1.8f}".format(axial)) + for turb in self.farm.turbines: + turb.axial = axial + turb.maxial.assign(axial) + adj_stop = time.time() + self.fprint("Axial Induction Adjusted: {:1.2f} s".format(adj_stop-adj_start),special="footer") def ChangeWindAngle(self,inflow_angle): """ @@ -352,12 +361,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 +376,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 +492,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") @@ -434,13 +525,15 @@ 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.000001) + eps=Constant(self.stability_eps) + eps.rename("eps","eps") nu = self.viscosity self.fprint("Viscosity: {:1.2e}".format(float(self.viscosity))) self.fprint("Max Mixing Length: {:1.2e}".format(float(self.lmax))) + rot_mat = as_tensor([[cos(inflow_angle),-sin(inflow_angle)],[sin(inflow_angle), cos(inflow_angle)]]) ### Create the test/trial/functions ### self.v,self.q = TestFunctions(self.fs.W) @@ -469,13 +562,14 @@ def ComputeFunctional(self,inflow_angle): 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 += - inner(dot(rot_mat,self.bd.bf_bk),self.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 in the Stabilizing term ### + if abs(float(eps)) >= 1e-14: + self.fprint("Using Stabilization Term") + stab = - eps*inner(grad(self.q), grad(self.p_k))*dx - eps*inner(grad(self.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: diff --git a/windse/RefinementManager.py b/windse/RefinementManager.py index 5384c83..ce9433b 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 c76dff9..04479b1 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,21 +163,30 @@ 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 self.problem.dom.mesh.coordinates()[:]=self.problem.dom.mesh.coordinates()[:]*self.problem.dom.xscale + def ChangeAxial(self,axial): + """ + This function recomputes all necessary components for a new axial induction + + Args: + axial (float): The new axial induction factor + """ + self.problem.ChangeAxial(axial) + def ChangeWindSpeed(self,inflow_speed): """ This function recomputes all necessary components for a new wind direction @@ -228,6 +238,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 @@ -458,10 +469,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) @@ -2066,6 +2077,7 @@ def __init__(self,problem): def Solve(self): for i, theta in enumerate(self.angles): + tic = time.time() self.fprint("Performing Solve {:d} of {:d}".format(i+1,len(self.angles)),special="header") self.fprint("Wind Angle: "+repr(math_to_meteor(theta))) if i > 0 or not near(theta,self.problem.dom.inflow_angle): @@ -2073,7 +2085,13 @@ def Solve(self): self.ChangeWindAngle(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") + + + toc = time.time() + + self.fprint("Finished Solve {:d} of {:d}: {:0.2f} s".format(i+1,len(self.angles),toc-tic),special="footer") + + class TimeSeriesSolver(SteadySolver): """ @@ -2091,29 +2109,47 @@ def __init__(self,problem): raise ValueError("Cannot use a Multi-Angle Solver with an "+self.params["domain"]["type"]+" domain.") self.orignal_solve = super(TimeSeriesSolver, self).Solve - raw_data = np.loadtxt(self.velocity_path,comments="#") + raw_data = np.genfromtxt(self.velocity_path,comments="#", delimiter=", ", skip_header=1) self.times = raw_data[:,0] self.speeds = raw_data[:,1] self.angles = raw_data[:,2] + self.axials = raw_data[:,4] self.num_solve = len(self.speeds) def Solve(self): + import time + for i in range(self.num_solve): + + tic = time.time() + time = self.times[i] theta = self.angles[i] speed = self.speeds[i] + axial = self.axials[i] self.fprint("Performing Solve {:d} of {:d}".format(i+1,len(self.angles)),special="header") self.fprint("Time: "+repr(time)) self.fprint("Wind Angle: "+repr(theta)) self.fprint("Wind Speed: "+repr(speed)) + self.fprint("Axial Val: "+repr(axial)) + + # TODO: need to handle the case where axial induction is zero + + # update values + self.ChangeAxial(axial) if i > 0 or not near(speed,self.problem.bd.HH_vel): self.ChangeWindSpeed(speed) if i > 0 or not near(theta,self.problem.dom.inflow_angle): self.ChangeWindAngle(theta) self.iter_val = time + + # solve self.orignal_solve() - self.fprint("Finished Solve {:d} of {:d}".format(i+1,len(self.angles)),special="footer") + toc = time.time() + + + self.fprint("Finished Solve {:d} of {:d}: {:0.2f} s".format(i+1,len(self.angles),toc-tic),special="footer") diff --git a/windse/__init__.py b/windse/__init__.py index 785bc91..72435a5 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/helper_functions.py b/windse/helper_functions.py index 7a077de..0010e9a 100644 --- a/windse/helper_functions.py +++ b/windse/helper_functions.py @@ -48,7 +48,7 @@ def math_to_meteor(angle): return new_angle - return 270.0-np.degrees(angle) + # return 270.0-np.degrees(angle) def ufl_eval(form, print_statement=None): ''' diff --git a/windse/input_schema.yaml b/windse/input_schema.yaml index bb0356d..155e7dc 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" @@ -855,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: + - "array" + description: "location in rotor diameters of the neighboring farm" # ================================================================ problem: # additionalProperties: False diff --git a/windse/wind_farm_types/GenericWindFarm.py b/windse/wind_farm_types/GenericWindFarm.py index 3b41287..4041c3e 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 @@ -420,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): """ diff --git a/windse/wind_farm_types/GridWindFarm.py b/windse/wind_farm_types/GridWindFarm.py index 4a954f6..8b1b07c 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/TrapWindFarm.py b/windse/wind_farm_types/TrapWindFarm.py new file mode 100644 index 0000000..e170aba --- /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 diff --git a/windse/wind_farm_types/__init__.py b/windse/wind_farm_types/__init__.py index a508f8c..c9b5f98 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.py b/windse_driver/driver.py index c492bd6..1b73a65 100644 --- a/windse_driver/driver.py +++ b/windse_driver/driver.py @@ -54,6 +54,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 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() ### Clean up other module references ### diff --git a/windse_driver/driver_functions.py b/windse_driver/driver_functions.py index cd66da3..acedb2e 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}