diff --git a/.github/workflows/CI_every_commit.yml b/.github/workflows/CI_every_commit.yml index 9e1ab959..4b993437 100644 --- a/.github/workflows/CI_every_commit.yml +++ b/.github/workflows/CI_every_commit.yml @@ -28,6 +28,7 @@ jobs: 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: | + conda install -n base mamba bash ./install.sh test-environment conda activate test-environment conda list diff --git a/demo/documented/Yaml_Examples/0-wind_farm_2D.yaml b/demo/documented/Yaml_Examples/0-wind_farm_2D.yaml index ee7a8d92..af7f84c6 100755 --- a/demo/documented/Yaml_Examples/0-wind_farm_2D.yaml +++ b/demo/documented/Yaml_Examples/0-wind_farm_2D.yaml @@ -60,11 +60,10 @@ domain: # ########################### Circle Domain ######################### type: circle - mesh_type: mshr # squircular, elliptic, stretch + mesh_type: gmsh # squircular, elliptic, stretch radius: 1200 # x-range of the domain | m center: [0.0, 0.0] # y-range of the domain | m nt: 62 # segments around circle| - - res: 50 # resolution for mshr | - # ################################################################### 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 a4d59c85..3b52b4f4 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 @@ -57,7 +57,7 @@ domain: ########################### Circle Domain ######################### type: circle - mesh_type: mshr # squircular, elliptic, stretch + mesh_type: gmsh # squircular, elliptic, stretch radius: 1500 # x-range of the domain | m center: [0.0, 0.0] # y-range of the domain | m nt: 50 # segments around circle| - 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 93a65a5a..8f51ac8e 100644 --- a/demo/documented/Yaml_Examples/7-wind_farm_2D_OM.yaml +++ b/demo/documented/Yaml_Examples/7-wind_farm_2D_OM.yaml @@ -18,8 +18,8 @@ wind_farm: type: grid # | grid_rows: 1 # Number of rows | - grid_cols: 2 # Number of columns | - - ex_x: [-882, 882] # x-extent of the farm | m - ex_y: [-882, 882] # y-extent of the farm | m + ex_x: [-720, 720] # x-extent of the farm | m + ex_y: [-720, 720] # y-extent of the farm | m jitter: 0 # Randomly perturb turbines| m seed: 8675309 # random seed for repeats | - ################################################################## @@ -57,11 +57,10 @@ domain: ########################### Circle Domain ######################### type: circle - mesh_type: mshr # squircular, elliptic, stretch + mesh_type: gmsh # squircular, elliptic, stretch radius: 1500 # x-range of the domain | m center: [0.0, 0.0] # y-range of the domain | m nt: 20 # segments around circle| - - res: 10 # resolution for mshr | - ################################################################### diff --git a/environment.yaml b/environment.yaml index 02730390..8fc0ec58 100644 --- a/environment.yaml +++ b/environment.yaml @@ -3,15 +3,17 @@ channels: - conda-forge - defaults dependencies: - # - python=3.7 + - python=3.9 #3.10 causes CI to hang in parallel + - mpich + # - mpich=4.1.2=hd33e60e_101 + # - mpich=4.0.3 - coveralls - dolfin-adjoint=2019.1.0 - - fenics=2019.1.0=py39*_26 + - fenics - hdf5 - jsonschema - matplotlib - memory_profiler - - mshr - openmdao - pandas - pyoptsparse @@ -19,7 +21,7 @@ dependencies: - pytest-cov - pytest-mpi - pyyaml - - scipy=1.7.0 + - scipy - slepc - sphinx<7.0.0 - pip @@ -34,4 +36,6 @@ dependencies: - networkx - pulp - fatpack + - gmsh + - meshio - -e . diff --git a/install.sh b/install.sh index 6d8594bd..1b2457f2 100644 --- a/install.sh +++ b/install.sh @@ -1,23 +1,15 @@ #!/bin/bash -### Run this to use conda in the script -source $(conda info --base)/etc/profile.d/conda.sh - -# ### Create the Environment -# conda create -y --name $1 -# conda activate $1 +### Create some names +conda_base_path=$(conda info --base) +env_name=$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 +### Run this to use conda in the script +source ${conda_base_path}/etc/profile.d/conda.sh -# ### Install editible version of WindSE -# pip install -e . +### Create the Environment with mamba +conda create --yes --name ${env_name} --channel conda-forge python=3.9 mamba +conda activate ${env_name} -conda env create --name ${1} --file environment.yaml +### Install dependencies from environment +mamba env update --prefix ${conda_base_path}/envs/${env_name} --file environment.yaml --prune diff --git a/tests/9-Regression/2D_Rectangle_Steady_gmsh.yaml b/tests/9-Regression/2D_Rectangle_Steady_gmsh.yaml new file mode 100644 index 00000000..3652e39d --- /dev/null +++ b/tests/9-Regression/2D_Rectangle_Steady_gmsh.yaml @@ -0,0 +1,81 @@ +# 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 + ################################################################### + +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: rectangle # | + x_range: [-1200, 1200] # x-range of the domain | m + y_range: [-600, 600] # y-range of the domain | m + nx: 128 # Number of x-nodes | - + ny: 64 # Number of y-nodes | - + mesh_type: gmsh + ################################################################### + +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: taylor_hood + +boundary_conditions: + vel_profile: uniform + 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: ["north","south"] + ################################################################## + +problem: + type: taylor_hood + viscosity: 5 + lmax: 50 + +solver: + # nonlinear_solver: newton + # newton_relaxation: 0.9 + type: steady + save_power: true + +optimization: + control_types: [layout,yaw] + layout_bounds: [[-720, 720],[-450, 450]] + gradient: True diff --git a/tests/9-Regression/3D_Box_Steady_gmsh.yaml b/tests/9-Regression/3D_Box_Steady_gmsh.yaml new file mode 100644 index 00000000..210c39f4 --- /dev/null +++ b/tests/9-Regression/3D_Box_Steady_gmsh.yaml @@ -0,0 +1,82 @@ +# 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 + ################################################################### + +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 | - + mesh_type: gmsh + ################################################################### + +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/2D_Rectangle_Steady_gmsh_truth.yaml b/tests/9-Regression/Truth_Data/2D_Rectangle_Steady_gmsh_truth.yaml new file mode 100644 index 00000000..b5f48e15 --- /dev/null +++ b/tests/9-Regression/Truth_Data/2D_Rectangle_Steady_gmsh_truth.yaml @@ -0,0 +1,108 @@ +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: 2 + x_min: -1200.0 + x_max: 1200.0 + y_min: -600.0 + y_max: 600.0 + volume: 2879999.9999999944 + ID_east: 1 + SA_east: 1199.9999999999993 + ID_north: 2 + SA_north: 2399.999999999999 + ID_west: 3 + SA_west: 1199.9999999999998 + ID_south: 4 + SA_south: 2399.999999999999 +FunctionSpaceManager: + velocity_dofs: 199784 + pressure_dofs: 25045 + total_dofs: 224829 +BoundaryManager: + min_x: 8.0 + max_x: 8.0 + avg_x: 8.0 + min_y: 0.0 + max_y: 0.0 + avg_y: 0.0 + min_p: 0.0 + max_p: 0.0 + avg_p: 0.0 + min_initial_values: 0.0 + max_initial_values: 8.0 + avg_initial_values: 3.5544169124089864 + num_bc: 7 +ProblemManager: + int_nu_T: 3.6195476564638705e-12 + int_tf_x: -0.017779690602830505 + int_tf_y: -0.006469952285041019 +SolverManager: + min_x_vel: 3.1513290589444534 + max_x_vel: 10.314223723630104 + avg_x_vel: 7.999999999999985 + min_y_vel: -1.5355717615335474 + max_y_vel: 1.7312382083757618 + avg_y_vel: -0.012640362247818808 +OptimizationManager: + n_controls: 24 + obj_value0: 0.09657585867884727 + 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: -8.33214537992983e-06 + grad_y_0: -2.0118478496424273e-05 + grad_x_1: -5.802770168646045e-07 + grad_y_1: -6.845843664967294e-07 + grad_x_2: 2.423778555536142e-06 + grad_y_2: 1.1422838544789697e-05 + grad_x_3: 8.578153855238626e-06 + grad_y_3: 1.6946275983151124e-05 + grad_x_4: -1.0314760814733896e-05 + grad_y_4: 2.8272220694531887e-05 + grad_x_5: -8.751615197887761e-07 + grad_y_5: -8.281878532430593e-06 + grad_x_6: 1.6008364572576094e-06 + grad_y_6: -1.2717509235006934e-05 + grad_x_7: 6.279591155513218e-06 + grad_y_7: -1.6236312575846403e-05 + grad_yaw_0: -0.0072120688648149335 + grad_yaw_1: -0.004558722084702321 + grad_yaw_2: -0.004826849442760917 + grad_yaw_3: -0.006876464533418879 + grad_yaw_4: -0.007472719447171188 + grad_yaw_5: -0.004770395674441498 + grad_yaw_6: -0.004642931548333499 + grad_yaw_7: -0.005883549707396079 diff --git a/tests/9-Regression/Truth_Data/3D_Box_Steady_gmsh_truth.yaml b/tests/9-Regression/Truth_Data/3D_Box_Steady_gmsh_truth.yaml new file mode 100644 index 00000000..018cd4a9 --- /dev/null +++ b/tests/9-Regression/Truth_Data/3D_Box_Steady_gmsh_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.9999719 + ID_east: 1 + SA_east: 767952.0000000016 + ID_north: 2 + SA_north: 1535904.000000002 + ID_west: 3 + SA_west: 767952.0000000014 + ID_south: 4 + SA_south: 1535904.0000000035 + ID_bottom: 5 + SA_bottom: 2879999.9999999953 + ID_top: 6 + SA_top: 2880000.000000008 +FunctionSpaceManager: + velocity_dofs: 13668 + pressure_dofs: 4556 + total_dofs: 18224 +BoundaryManager: + min_x: 3.4975322130308117 + max_x: 9.14464606512413 + avg_x: 7.757844631064229 + 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.9394611577660588 + num_bc: 11 +ProblemManager: + int_nu_T: 2.6163033631208665 + int_tf_x: -0.0028305417152973003 + int_tf_y: -0.0010300218517793624 + int_tf_z: 0.0 +SolverManager: + min_x_vel: 0.0 + max_x_vel: 10.040397540792812 + avg_x_vel: 8.170579698847535 + min_y_vel: -0.9531449010214702 + max_y_vel: 0.6420538617788053 + avg_y_vel: 0.012192747697076073 + min_z_vel: -0.48671715986568975 + max_z_vel: 0.701965608955228 + avg_z_vel: 0.012192747697076073 +OptimizationManager: + n_controls: 24 + obj_value0: 6.323599124013655 + 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.004755930209050213 + grad_y_0: 0.0017881918659346323 + grad_x_1: 0.0003659875531067487 + grad_y_1: 0.0005655959904036023 + grad_x_2: 0.0004083746380206747 + grad_y_2: 0.00044901198168019853 + grad_x_3: 0.00039986116732627225 + grad_y_3: 0.0002979237791272051 + grad_x_4: 0.00023192273756160627 + grad_y_4: -0.0022518314669599535 + grad_x_5: -0.0027744685842574822 + grad_y_5: -0.001288273590947986 + grad_x_6: 0.0007418418324504838 + grad_y_6: 0.0018555755795040614 + grad_x_7: 0.0007795194897304387 + grad_y_7: 0.0017116373798768318 + grad_yaw_0: -0.6239086775156708 + grad_yaw_1: -0.4108354430166945 + grad_yaw_2: -0.4340509099311657 + grad_yaw_3: -0.44622355672999786 + grad_yaw_4: -0.6765518143158551 + grad_yaw_5: -0.32162601141179337 + grad_yaw_6: -0.37027251882928364 + grad_yaw_7: -0.5014367696212912 diff --git a/windse/DomainManager.py b/windse/DomainManager.py index 6279b5bc..e14e538a 100644 --- a/windse/DomainManager.py +++ b/windse/DomainManager.py @@ -1,5 +1,5 @@ """ -The DomainManager submodule contains the various classes used for +The DomainManager submodule contains the various classes used for creating different types of domains """ @@ -17,8 +17,6 @@ ### This checks if we are just doing documentation ### if not main_file in ["sphinx-build", "__main__.py"]: from dolfin import * - from mshr import * - #from windse.mshr_fake import * import copy import time import warnings @@ -35,7 +33,7 @@ ### Check if we need dolfin_adjoint ### if windse_parameters.dolfin_adjoint: from dolfin_adjoint import * - + ### This import improves the plotter functionality on Mac ### if platform == 'darwin': import matplotlib @@ -55,10 +53,10 @@ class SinglePeriodicBoundary(SubDomain): def __init__(self,a,b,axis=0): super(SinglePeriodicBoundary, self).__init__() self.axis = axis # 0 for periodic x and 1 for periodic y - self.target = a # location of the primary boundary - self.offset = b-a # the offset from the primary to the secondary + self.target = a # location of the primary boundary + self.offset = b-a # the offset from the primary to the secondary - # returns true primary boundary + # returns true primary boundary def inside(self, x, on_boundary): return near(x[self.axis],self.target) and on_boundary @@ -71,13 +69,13 @@ class DoublePeriodicBoundary(SubDomain): def __init__(self,x_range,y_range): super(DoublePeriodicBoundary, self).__init__() - self.x_range = x_range # location of the primary boundary + self.x_range = x_range # location of the primary boundary self.y_range = y_range # location of the primary boundary - self.x_offset = x_range[1]-x_range[0] # the offset from the primary to the secondary - self.y_offset = y_range[1]-y_range[0] # the offset from the primary to the secondary + self.x_offset = x_range[1]-x_range[0] # the offset from the primary to the secondary + self.y_offset = y_range[1]-y_range[0] # the offset from the primary to the secondary self.eps = 1e-4 - # returns true primary boundary + # returns true primary boundary def inside(self, x, on_boundary): on_either = near(x[0], self.x_range[0],self.eps) or near(x[1], self.y_range[0],self.eps) on_corner1 = near(x[0], self.x_range[1],self.eps) and near(x[1], self.y_range[0],self.eps) @@ -232,7 +230,7 @@ def DebugOutput(self): def Plot(self): """ - This function plots the domain using matplotlib and saves the + This function plots the domain using matplotlib and saves the output to output/.../plots/mesh.pdf """ @@ -395,7 +393,7 @@ def CylinderRefine(self,center,radius,height=0,expand_factor=1): self.Refine(cell_f) refine_stop = time.time() - self.fprint("Mesh Refinement Finished: {:1.2f} s".format(refine_stop-refine_start),special="footer") + self.fprint("Mesh Refinement Finished: {:1.2f} s".format(refine_stop-refine_start),special="footer") def StreamRefine(self,center,radius,length,theta=0,pivot_offset=0,expand_factor=1): refine_start = time.time() @@ -443,7 +441,7 @@ def StreamRefine(self,center,radius,length,theta=0,pivot_offset=0,expand_factor= self.Refine(cell_f) refine_stop = time.time() - self.fprint("Mesh Refinement Finished: {:1.2f} s".format(refine_stop-refine_start),special="footer") + self.fprint("Mesh Refinement Finished: {:1.2f} s".format(refine_stop-refine_start),special="footer") def Refine(self, cellmarkers=None): @@ -503,15 +501,15 @@ def Refine(self, cellmarkers=None): # num (int): the number of times to refine # :Keyword Arguments: - # * **region** (*list*): + # * **region** (*list*): # for square region use: [[xmin,xmax],[ymin,ymax],[zmin,zmax]] # for circle region use: [[radius],[c_x,c_y],[zmin,zmax]] - # * **region_type** (*str*): Either "circle" or "square" - # * **cell_markers** (*:meth:dolfin.mesh.MeshFunction*): A cell function marking which cells to refine + # * **region_type** (*str*): Either "circle" or "square" + # * **cell_markers** (*:meth:dolfin.mesh.MeshFunction*): A cell function marking which cells to refine # """ # refine_start = time.time() - # ### Print some useful stats + # ### Print some useful stats # self.fprint("Starting Mesh Refinement",special="header") # if cell_markers is not None: # self.fprint("Region Type: {0}".format("cell_markers")) @@ -585,7 +583,7 @@ def Refine(self, cellmarkers=None): # cells_marked += 1 # self.fprint("Cells Marked for Refinement: {:d}".format(cells_marked)) - + # ### If neither a region or cell markers were provided, Refine everwhere ### # else: # cell_f = MeshFunction('bool', self.mesh, self.mesh.geometry().dim(),True) @@ -604,15 +602,15 @@ def Refine(self, cellmarkers=None): # if num>1: # step_stop = time.time() # self.fprint("Step {:d} of {:d} Finished: {:1.2f} s".format(i+1,num,step_stop-step_start), special="footer") - + # refine_stop = time.time() # self.fprint("Mesh Refinement Finished: {:1.2f} s".format(refine_stop-refine_start),special="footer") def WarpSplit(self,h,s): """ - This function warps the mesh to shift more cells towards the ground. - is achieved by spliting the domain in two and moving the cells so + This function warps the mesh to shift more cells towards the ground. + is achieved by spliting the domain in two and moving the cells so that a percentage of them are below the split. Args: @@ -659,7 +657,7 @@ def WarpSplit(self,h,s): def WarpSmooth(self,s): """ - This function warps the mesh to shift more cells towards the ground. + This function warps the mesh to shift more cells towards the ground. The cells are shifted based on the function: .. math:: @@ -710,7 +708,7 @@ def transform(x,y,z,z0,z1): self.bmesh.coordinates()[:,1]=y_hd self.bmesh.coordinates()[:,2]=z_hd self.bmesh.bounding_box_tree().build(self.bmesh) - + self.fprint("Moving Mesh to New Boundary Using ALE") ### The way dolfin_adjoint overwrites ALE.move is destructive so we have to use a work around if self.params.dolfin_adjoint: @@ -897,12 +895,12 @@ def RecomputeBoundaryMarkers(self,inflow_angle): class BoxDomain(GenericDomain): """ A box domain is simply a 3D rectangular prism. This box is defined - by 6 parameters in the param.yaml file. + by 6 parameters in the param.yaml file. Example: In the yaml file define:: - domain: + domain: # # Description | Units x_range: [-2500, 2500] # x-range of the domain | m y_range: [-2500, 2500] # y-range of the domain | m @@ -911,7 +909,7 @@ class BoxDomain(GenericDomain): ny: 10 # Number of y-nodes | - nz: 2 # Number of z-nodes | - - This will produce a box with corner points (-2500,-2500,0.04) + This will produce a box with corner points (-2500,-2500,0.04) to (2500,2500,630). The mesh will have *nx* nodes in the *x*-direction, *ny* in the *y*-direction, and *nz* in the *z*-direction. """ @@ -937,9 +935,129 @@ def __init__(self): mesh_start = time.time() self.fprint("") self.fprint("Generating Mesh") - start = Point(self.x_range[0], self.y_range[0], self.z_range[0]) - stop = Point(self.x_range[1], self.y_range[1], self.z_range[1]) - self.mesh = BoxMesh(start, stop, self.nx, self.ny, self.nz) + + if self.mesh_type == "gmsh": + + if (self.params.rank == 0): + + # only dependencies for the gmsh case + import gmsh + import meshio + import tempfile + + # start up gmsh model + gmsh.initialize() + gmsh.model.add("boxmesh") + + # true extents + Lx_true = self.x_range[1] - self.x_range[0] + Ly_true = self.y_range[1] - self.y_range[0] + Lz_true = self.z_range[1] - self.z_range[0] + + # either do stretching to preserve + stretch_anisotropic = True # TODO: turn into an input param + if stretch_anisotropic: + Lx_prime = 1.0 + Ly_prime = self.ny/self.nx + Lz_prime = self.nz/self.nx + else: + Lx_prime = Lx_true + Ly_prime = Ly_true + Lz_prime = Lz_true + + # create points that constitute box's lower surface + pt_bx000 = gmsh.model.geo.addPoint( + self.x_range[0]/Lx_true*Lx_prime, + self.y_range[0]/Ly_true*Ly_prime, + self.z_range[0]/Lz_true*Lz_prime, + ) + pt_bx100 = gmsh.model.geo.addPoint( + self.x_range[1]/Lx_true*Lx_prime, + self.y_range[0]/Ly_true*Ly_prime, + self.z_range[0]/Lz_true*Lz_prime, + ) + pt_bx010 = gmsh.model.geo.addPoint( + self.x_range[0]/Lx_true*Lx_prime, + self.y_range[1]/Ly_true*Ly_prime, + self.z_range[0]/Lz_true*Lz_prime, + ) + pt_bx110 = gmsh.model.geo.addPoint( + self.x_range[1]/Lx_true*Lx_prime, + self.y_range[1]/Ly_true*Ly_prime, + self.z_range[0]/Lz_true*Lz_prime, + ) + + # create lines of box's lower surface + ln_bx000 = gmsh.model.geo.addLine(pt_bx000, pt_bx100) + ln_bx001 = gmsh.model.geo.addLine(pt_bx100, pt_bx110) + ln_bx010 = gmsh.model.geo.addLine(pt_bx010, pt_bx110) + ln_bx011 = gmsh.model.geo.addLine(pt_bx000, pt_bx010) + + # create curve loop for lower surface + cl_bx0 = gmsh.model.geo.addCurveLoop( + [ln_bx000, ln_bx001, -ln_bx010, -ln_bx011] + ) + # create plane surface + ps_bx0 = gmsh.model.geo.addPlaneSurface([cl_bx0]) + + bx_simple = gmsh.model.geo.extrude( + [(2, ps_bx0)], + 0.0, + 0.0, + (self.z_range[1] - self.z_range[0])/Lz_true*Lz_prime + ) + + # synchronize the geometry + gmsh.model.geo.synchronize() + + # characteristic length based on regular tetrahedron volume formula + lc_mean = ( + 6*np.sqrt(2) + *(Lx_prime/self.nx) + *(Ly_prime/self.ny) + *(Lz_prime/self.nz) + )**(1/3) + # set the mesh size to the characteristic length + gmsh.model.mesh.setSize(gmsh.model.getEntities(0), lc_mean) + + # generate and output + gmsh.model.mesh.generate(3) + + # create a temporary directory to hold mesh IO files for conversion + dir_for_meshes = tempfile.TemporaryDirectory() + + # write and finalize gmsh + gmsh.write(os.path.join(dir_for_meshes.name, "dummy.msh")) + gmsh.finalize() # don't need it anymore! + + # use meshio to convert from gmsh to dolfin xml file + mesh_meshio = meshio.read(os.path.join(dir_for_meshes.name, "dummy.msh")) + meshio.write( + os.path.join(dir_for_meshes.name, "dummy.xml"), + mesh_meshio, + file_format="dolfin-xml", + ) + + # broadcast the temporary directory for shared use + dirname_for_meshes = self.params.comm.bcast(dir_for_meshes.name, root=0) + + # load the xml into windse + self.mesh = Mesh(os.path.join(dirname_for_meshes, "dummy.xml")) + + # delete the temporary directory + if self.params.rank == 0: + dir_for_meshes.cleanup() + + # stretch the coordinates back to the physical dimensions (now with the + # specified aspect ratios) + self.mesh.coordinates()[:,0] = self.mesh.coordinates()[:,0]*Lx_true/Lx_prime + self.mesh.coordinates()[:,1] = self.mesh.coordinates()[:,1]*Ly_true/Ly_prime + self.mesh.coordinates()[:,2] = self.mesh.coordinates()[:,2]*Lz_true/Lz_prime + + else: + start = Point(self.x_range[0], self.y_range[0], self.z_range[0]) + stop = Point(self.x_range[1], self.y_range[1], self.z_range[1]) + self.mesh = BoxMesh(start, stop, self.nx, self.ny, self.nz) # box = Box(start,stop) # self.mesh = generate_mesh(box,self.nx) @@ -992,7 +1110,7 @@ def RecomputeBoundaryMarkers(self,inflow_angle): north_id = self.boundary_names["north"] west_id = self.boundary_names["west"] south_id = self.boundary_names["south"] - + ### Set up the baseline (angle=0) order ### cardinal_ids = [east_id,north_id,west_id,south_id] diagonal_ids = [outflow_id,outflow_id,inflow_id,inflow_id] @@ -1010,7 +1128,7 @@ def RecomputeBoundaryMarkers(self,inflow_angle): ### Get a list of all wall facets ### wall_facets = [] for i in cardinal_ids: - wall_facets += self.boundary_markers.where_equal(i) + wall_facets += self.boundary_markers.where_equal(i) ### Iterate through facets remarking as prescribed by new_order ### for facet_id in wall_facets: @@ -1032,22 +1150,22 @@ def RecomputeBoundaryMarkers(self,inflow_angle): class CylinderDomain(GenericDomain): """ - A cylinder domain is a cylinder that is centered a c0 and has radius r. - This domain is defined by 6 parameters in the param.yaml file. The + A cylinder domain is a cylinder that is centered a c0 and has radius r. + This domain is defined by 6 parameters in the param.yaml file. The center of the cylinder is assumed to be the z-axis. Example: In the yaml file define:: - domain: + domain: # # Description | Units z_range: [0.04, 630] # z-range of the domain | m radius: 2500 # radius of base circle | m nt: 100 # Number of radial nodes| - nz: 10 # Number of z nodes | - - This will produce a upright cylinder centered at (0.0,0.0) with a - radius of 2500 m and extends from z=0.04 to 630 m. The mesh will + This will produce a upright cylinder centered at (0.0,0.0) with a + radius of 2500 m and extends from z=0.04 to 630 m. The mesh will have *nx* nodes in the *x*-direction, *ny* in the *y*-direction, and *nz* in the *z*-direction. """ @@ -1077,22 +1195,105 @@ def __init__(self): mesh_start = time.time() self.fprint("") if self.mesh_type == "mshr": - self.fprint("Generating Mesh Using mshr") + raise NotImplementedError("Mshr is no longer supported, use gmsh instead.") + # self.fprint("Generating Mesh Using mshr") + + # ### Create Mesh ### + # # mshr_circle = Circle(Point(self.center[0],self.center[1]), self.radius, self.nt) + # # mshr_domain = Extrude2D(mshr_circle,self.z_range[1]-self.z_range[0]) + # top = Point(self.center[0],self.center[1],self.z_range[1]) + # bottom = Point(self.center[0],self.center[1],self.z_range[0]) + # mshr_domain = Cylinder(top,bottom,self.radius,self.radius,self.nt) + # self.mesh = generate_mesh(mshr_domain,self.res) + # # self.mesh = refine(self.mesh) + # # self.mesh = refine(self.mesh) + # # self.mesh = refine(self.mesh) + + # # z = self.mesh.coordinates()[:,2]#+self.z_range[0] + # # self.mesh.coordinates()[:,2] = z + # # self.mesh.bounding_box_tree().build(self.mesh) + + elif self.mesh_type == "gmsh": + self.fprint("Generating Mesh Using gmsh") + + if (self.params.rank == 0): + + # only dependencies for the gmsh case + import gmsh + import meshio + import tempfile + + # start up gmsh model + gmsh.initialize() + gmsh.model.add("circmesh") + factory = gmsh.model.occ + + # true extents + R_true = self.radius + Lz_true = self.z_range[1] - self.z_range[0] + + stretch_anisotropic = True # TODO: turn into an input param + if stretch_anisotropic: + R_prime = 1.0 + Lz_prime = 2*np.pi*self.nz/self.nt + + lc_mean = Lz_prime/self.nz + else: + raise NotImplementedError("logic for ignoring scaling not implemented. -cfrontin") - ### Create Mesh ### - # mshr_circle = Circle(Point(self.center[0],self.center[1]), self.radius, self.nt) - # mshr_domain = Extrude2D(mshr_circle,self.z_range[1]-self.z_range[0]) - top = Point(self.center[0],self.center[1],self.z_range[1]) - bottom = Point(self.center[0],self.center[1],self.z_range[0]) - mshr_domain = Cylinder(top,bottom,self.radius,self.radius,self.nt) - self.mesh = generate_mesh(mshr_domain,self.res) - # self.mesh = refine(self.mesh) - # self.mesh = refine(self.mesh) - # self.mesh = refine(self.mesh) - - # z = self.mesh.coordinates()[:,2]#+self.z_range[0] - # self.mesh.coordinates()[:,2] = z - # self.mesh.bounding_box_tree().build(self.mesh) + # set up the geometry + x_circ = 0.0 + y_circ = 0.0 + z_circ = 0.0 # self.z_range[0] + + # create a circle in the horizon plane + ln_circ = factory.addCircle(x_circ, y_circ, z_circ, R_prime) + cl_circ = factory.addCurveLoop([ln_circ]) + ps_circ = factory.addPlaneSurface([cl_circ]) + + geom_cyl = factory.extrude([(2, ps_circ)], 0.0, 0.0, Lz_prime) + + # synchronize the geometry + factory.synchronize() + + # set the mesh size + gmsh.model.mesh.setSize(gmsh.model.getEntities(0), lc_mean) + + # generate and output + gmsh.model.mesh.generate(3) + # gmsh.write("dummy.msh") # DEBUG!!!!! + + # create a temporary directory to hold mesh IO files for conversion + dir_for_meshes = tempfile.TemporaryDirectory() + + # write and finalize gmsh + gmsh.write(os.path.join(dir_for_meshes.name, "dummy.msh")) + gmsh.finalize() # don't need it anymore! + + # use meshio to convert from gmsh to dolfin xml file + mesh_meshio = meshio.read(os.path.join(dir_for_meshes.name, "dummy.msh")) + + meshio.write( + os.path.join(dir_for_meshes.name, "dummy.xml"), + mesh_meshio, + file_format="dolfin-xml", + ) + + # broadcast the temporary directory for shared use + dirname_for_meshes = self.params.comm.bcast(dir_for_meshes.name, root=0) + + # load the xml into windse + self.mesh = Mesh(os.path.join(dirname_for_meshes, "dummy.xml")) + + # delete the temporary directory + if self.params.rank == 0: + dir_for_meshes.cleanup() + + # stretch the coordinates back to the physical dimensions (now with the + # specified aspect ratios) + self.mesh.coordinates()[:,0] = self.mesh.coordinates()[:,0]*R_true/R_prime + self.center[0] + self.mesh.coordinates()[:,1] = self.mesh.coordinates()[:,1]*R_true/R_prime + self.center[1] + self.mesh.coordinates()[:,2] = self.mesh.coordinates()[:,2]*Lz_true/Lz_prime + self.z_range[0] else: self.fprint("Generating Box Mesh") @@ -1107,7 +1308,7 @@ def __init__(self): x = self.mesh.coordinates()[:,0] y = self.mesh.coordinates()[:,1] z = self.mesh.coordinates()[:,2] - + self.fprint("Morphing Mesh") if self.mesh_type == "elliptic": x_hat, y_hat, z_hat = Elliptical_Grid(x, y, z, self.radius) @@ -1136,10 +1337,29 @@ def __init__(self): mark_start = time.time() self.fprint("") self.fprint("Marking Boundaries") - outflow = CompiledSubDomain("on_boundary", nx=nom_x, ny=nom_y, z0 = self.z_range[0], z1 = self.z_range[1]) - inflow = CompiledSubDomain("nx*(x[0]-c0)+ny*(x[1]-c1)<=0 && on_boundary", nx=nom_x, ny=nom_y, z0 = self.z_range[0], z1 = self.z_range[1], c0=self.center[0], c1=self.center[1]) - top = CompiledSubDomain("near(x[2], z1, tol) && on_boundary",z1 = self.z_range[1],tol = 1e-10) - bottom = CompiledSubDomain("near(x[2], z0, tol) && on_boundary",z0 = self.z_range[0],tol = 1e-10) + outflow = CompiledSubDomain( + "on_boundary", + nx=nom_x, ny=nom_y, + z0=self.z_range[0], z1=self.z_range[1], + ) + inflow = CompiledSubDomain( + "nx*(x[0]-c0)+ny*(x[1]-c1)<=0 && on_boundary", + nx=nom_x, ny=nom_y, + z0=self.z_range[0], z1=self.z_range[1], + c0=self.center[0], c1=self.center[1], + ) + top = CompiledSubDomain( + "near(x[2], z1, tol) && on_boundary", + z1=self.z_range[1], tol=1e-10 + ) + bottom = CompiledSubDomain( + "near(x[2], z0, tol) && on_boundary", + z0=self.z_range[0], tol=1e-10, + ) + # outflow = CompiledSubDomain("x[2] > z0 && x[2] < z1 && on_boundary", nx=nom_x, ny=nom_y, z0 = self.z_range[0], z1 = self.z_range[1]) + # inflow = CompiledSubDomain("x[2] > z0 && x[2] < z1 && nx*(x[0]-c0)+ny*(x[1]-c1)<=0 && on_boundary", nx=nom_x, ny=nom_y, z0 = self.z_range[0], z1 = self.z_range[1], c0=self.center[0], c1=self.center[1]) + # top = CompiledSubDomain("near(x[2], z1, tol) && on_boundary",z1 = self.z_range[1],tol = 1e-10) + # bottom = CompiledSubDomain("near(x[2], z0, tol) && on_boundary",z0 = self.z_range[0],tol = 1e-10) self.boundary_subdomains = [None,None,None,None,outflow,inflow,bottom,top] self.boundary_names = {"west":None,"east":None,"south":None,"north":None,"bottom":7,"top":8,"inflow":6,"outflow":5} self.boundary_types = {"inflow": ["inflow"], @@ -1234,12 +1454,81 @@ def __init__(self): mesh_start = time.time() self.fprint("") if self.mesh_type == "mshr": + raise NotImplementedError("Mshr is no longer supported, use gmsh instead.") + # self.fprint("Generating Mesh Using mshr") + + # ### Create Mesh ### + # mshr_circle = Circle(Point(self.center[0],self.center[1]), self.radius, self.nt) + # self.mesh = generate_mesh(mshr_circle,self.res) + + + elif self.mesh_type == "gmsh": + self.fprint("Generating Mesh Using gmsh") + + if (self.params.rank == 0): + + # only dependencies for the gmsh case + import gmsh + import meshio + import tempfile + + # start up gmsh model + gmsh.initialize() + gmsh.model.add("circmesh") + factory = gmsh.model.occ + + # use arc division to set characteristic length on the horizon + lc_mean = 2*np.pi*self.radius/self.nt + + # set up the geometry + x_circ = self.center[0] + y_circ = self.center[1] + z_circ = 0.0 + + # create the circle in the horizon plane + ln_circ = factory.addCircle(x_circ, y_circ, z_circ, self.radius) + cl_circ = factory.addCurveLoop([ln_circ]) + ps_circ = factory.addPlaneSurface([cl_circ]) + + # synchronize the geometry + factory.synchronize() + + # set the mesh size + gmsh.model.mesh.setSize(gmsh.model.getEntities(0), lc_mean) + + # generate and output + gmsh.model.mesh.generate(2) + + # create a temporary directory to hold mesh IO files for conversion + dir_for_meshes = tempfile.TemporaryDirectory() + + # write and finalize gmsh + gmsh.write(os.path.join(dir_for_meshes.name, "dummy.msh")) + gmsh.finalize() # don't need it anymore! - self.fprint("Generating Mesh Using mshr") + # use meshio to convert from gmsh to dolfin xml file + mesh_meshio = meshio.read(os.path.join(dir_for_meshes.name, "dummy.msh")) - ### Create Mesh ### - mshr_circle = Circle(Point(self.center[0],self.center[1]), self.radius, self.nt) - self.mesh = generate_mesh(mshr_circle,self.res) + # flatten to 2D + assert np.all(np.isclose(mesh_meshio.points[0,2], mesh_meshio.points[:,2])) + mesh_meshio.points = mesh_meshio.points[:,0:2] + + # write the dolfin file, and then finally re-load it into dolfin + meshio.write( + os.path.join(dir_for_meshes.name, "dummy.xml"), + mesh_meshio, + file_format="dolfin-xml", + ) + + # broadcast the temporary directory for shared use + dirname_for_meshes = self.params.comm.bcast(dir_for_meshes.name, root=0) + + # load the xml into windse + self.mesh = Mesh(os.path.join(dirname_for_meshes, "dummy.xml")) + + # delete the temporary directory + if self.params.rank == 0: + dir_for_meshes.cleanup() else: self.fprint("Generating Rectangle Mesh") @@ -1253,7 +1542,7 @@ def __init__(self): self.mesh = refine(self.mesh) x = self.mesh.coordinates()[:,0] y = self.mesh.coordinates()[:,1] - + self.fprint("Morphing Mesh") if self.mesh_type == "elliptic": x_hat, y_hat, z_hat = Elliptical_Grid(x, y, 0, self.radius) @@ -1352,19 +1641,19 @@ def RecomputeBoundaryMarkers(self,inflow_angle): class RectangleDomain(GenericDomain): """ A rectangle domain is simply a 2D rectangle. This mesh is defined - by 4 parameters in the param.yaml file. + by 4 parameters in the param.yaml file. Example: In the yaml file define:: - domain: + domain: # # Description | Units x_range: [-2500, 2500] # x-range of the domain | m y_range: [-2500, 2500] # y-range of the domain | m nx: 10 # Number of x-nodes | - ny: 10 # Number of y-nodes | - - This will produce a rectangle with corner points (-2500,-2500) + This will produce a rectangle with corner points (-2500,-2500) to (2500,2500). The mesh will have *nx* nodes in the *x*-direction, and *ny* in the *y*-direction. @@ -1388,9 +1677,118 @@ def __init__(self): mesh_start = time.time() self.fprint("") self.fprint("Generating Mesh") - start = Point(self.x_range[0], self.y_range[0]) - stop = Point(self.x_range[1], self.y_range[1]) - self.mesh = RectangleMesh(start, stop, self.nx, self.ny) + + if self.mesh_type == "gmsh": + + if (self.params.rank == 0): + + # only dependencies for the gmsh case + import gmsh + import meshio + import tempfile + + # start up gmsh model + gmsh.initialize() + gmsh.model.add("rectmesh") + + # true extents + Lx_true = self.x_range[1] - self.x_range[0] + Ly_true = self.y_range[1] - self.y_range[0] + + # either do stretching to preserve + stretch_anisotropic = True # TODO: turn into an input param + if stretch_anisotropic: + Lx_prime = 1.0 + Ly_prime = self.ny/self.nx + else: + Lx_prime = Lx_true + Ly_prime = Ly_true + + # create points that constitute box's lower surface + pt_bx000 = gmsh.model.geo.addPoint( + self.x_range[0]/Lx_true*Lx_prime, + self.y_range[0]/Ly_true*Ly_prime, + 0.0, + ) + pt_bx100 = gmsh.model.geo.addPoint( + self.x_range[1]/Lx_true*Lx_prime, + self.y_range[0]/Ly_true*Ly_prime, + 0.0, + ) + pt_bx010 = gmsh.model.geo.addPoint( + self.x_range[0]/Lx_true*Lx_prime, + self.y_range[1]/Ly_true*Ly_prime, + 0.0, + ) + pt_bx110 = gmsh.model.geo.addPoint( + self.x_range[1]/Lx_true*Lx_prime, + self.y_range[1]/Ly_true*Ly_prime, + 0.0, + ) + + # create lines of box's lower surface + ln_bx000 = gmsh.model.geo.addLine(pt_bx000, pt_bx100) + ln_bx001 = gmsh.model.geo.addLine(pt_bx100, pt_bx110) + ln_bx010 = gmsh.model.geo.addLine(pt_bx010, pt_bx110) + ln_bx011 = gmsh.model.geo.addLine(pt_bx000, pt_bx010) + + # create curve loop for lower surface + cl_bx0 = gmsh.model.geo.addCurveLoop( + [ln_bx000, ln_bx001, -ln_bx010, -ln_bx011] + ) + # create plane surface + ps_bx0 = gmsh.model.geo.addPlaneSurface([cl_bx0]) + + # synchronize the geometry + gmsh.model.geo.synchronize() + + # characteristic length based on regular tetrahedron volume formula + lc_mean = np.sqrt(4/np.sqrt(3)*(Lx_prime/self.nx)*(Ly_prime/self.ny)) + # set the mesh size to the characteristic length + gmsh.model.mesh.setSize(gmsh.model.getEntities(0), lc_mean) + + # generate and output + gmsh.model.mesh.generate(2) + + # create a temporary directory to hold mesh IO files for conversion + dir_for_meshes = tempfile.TemporaryDirectory() + + # write and finalize gmsh + gmsh.write(os.path.join(dir_for_meshes.name, "dummy.msh")) + gmsh.finalize() # don't need it anymore! + + # use meshio to convert from gmsh to dolfin xml file + mesh_meshio = meshio.read(os.path.join(dir_for_meshes.name, "dummy.msh")) + + # flatten to 2D + assert np.all(np.isclose(mesh_meshio.points[0,2], mesh_meshio.points[:,2])) + mesh_meshio.points = mesh_meshio.points[:,0:2] + + meshio.write( + os.path.join(dir_for_meshes.name, "dummy.xml"), + mesh_meshio, + file_format="dolfin-xml", + ) + + # broadcast the temporary directory for shared use + dirname_for_meshes = self.params.comm.bcast(dir_for_meshes.name, root=0) + + # load the xml into windse + self.mesh = Mesh(os.path.join(dirname_for_meshes, "dummy.xml")) + + # delete the temporary directory + if self.params.rank == 0: + dir_for_meshes.cleanup() + + # stretch the coordinates back to the physical dimensions (now with the + # specified aspect ratios) + self.mesh.coordinates()[:,0] = self.mesh.coordinates()[:,0]*Lx_true/Lx_prime + self.mesh.coordinates()[:,1] = self.mesh.coordinates()[:,1]*Ly_true/Ly_prime + else: + start = Point(self.x_range[0], self.y_range[0]) + stop = Point(self.x_range[1], self.y_range[1]) + self.mesh = RectangleMesh(start, stop, self.nx, self.ny) + self.bmesh = BoundaryMesh(self.mesh,"exterior") mesh_stop = time.time() self.fprint("Mesh Generated: {:1.2f} s".format(mesh_stop-mesh_start)) @@ -1431,7 +1829,7 @@ def RecomputeBoundaryMarkers(self,inflow_angle): north_id = self.boundary_names["north"] west_id = self.boundary_names["west"] south_id = self.boundary_names["south"] - + ### Set up the baseline (angle=0) order ### cardinal_ids = [east_id,north_id,west_id,south_id] diagonal_ids = [outflow_id,outflow_id,inflow_id,inflow_id] @@ -1449,7 +1847,7 @@ def RecomputeBoundaryMarkers(self,inflow_angle): ### Get a list of all wall facets ### wall_facets = [] for i in cardinal_ids: - wall_facets += self.boundary_markers.where_equal(i) + wall_facets += self.boundary_markers.where_equal(i) ### Iterate through facets remarking as prescribed by new_order ### for facet_id in wall_facets: @@ -1472,23 +1870,23 @@ def RecomputeBoundaryMarkers(self,inflow_angle): class ImportedDomain(GenericDomain): """ This class generates a domain from imported files. This mesh is defined - by 2 parameters in the param.yaml file. + by 2 parameters in the param.yaml file. Example: In the yaml file define:: - domain: + domain: path: "Mesh_data/" filetype: "xml.gz" The supported filetypes are "xml.gz" and "h5". For "xml.gz" 3 files are - required: + required: * mesh.xml.gz - this contains the mesh in a format dolfin can handle * boundaries.xml.gz - this contains the facet markers that define where the boundaries are - * topology.txt - this contains the data for the ground topology. + * topology.txt - this contains the data for the ground topology. It assumes that the coordinates are from a uniform mesh. - It contains three column: x, y, z. The x and y columns contain + It contains three column: x, y, z. The x and y columns contain just the unique values. The z column contains the ground values for every combination of x and y. The first row must be the number of points in the x and y direction. Here is an example for z=x+y/10:: diff --git a/windse/FunctionSpaceManager.py b/windse/FunctionSpaceManager.py index 66ee873a..ab78c00a 100644 --- a/windse/FunctionSpaceManager.py +++ b/windse/FunctionSpaceManager.py @@ -89,7 +89,7 @@ class LinearFunctionSpace(GenericFunctionSpace): def __init__(self,dom): super(LinearFunctionSpace, self).__init__(dom) - # trick the mesh to working? + # # trick the mesh to working? # dummy = MeshFunction('bool', self.mesh, self.mesh.geometry().dim(),False) # print("before:", self.mesh.num_entities_global(0)) # self.mesh = refine(self.mesh,dummy) diff --git a/windse/input_schema.yaml b/windse/input_schema.yaml index 18291dd0..d89c1c9f 100644 --- a/windse/input_schema.yaml +++ b/windse/input_schema.yaml @@ -18,7 +18,7 @@ properties: # - required_value properties: name: - type: + type: - "string" - "null" description: "The folder name for the current run." @@ -89,7 +89,6 @@ properties: - center - radius - nt - - res - interpolated - if: properties: @@ -101,7 +100,6 @@ properties: - center - radius - nt - - res - if: properties: type: @@ -122,9 +120,9 @@ properties: - terrain_path - interpolated - gaussian - properties: + properties: type: - type: + type: - "string" - "null" description: "The type of domain, options are ``rectangle``, ``box``, ``cylinder``, ``circle``, ``imported``, or ``interpolated``. See :ref:`imported_domain`" @@ -138,25 +136,25 @@ properties: - "interpolated" - Null path: - type: + type: - "string" - "null" description: "The path to all domain files if using imported domain and standard naming." default: Null mesh_path: - type: + type: - "string" - "null" description: "A specific path for the imported mesh file." default: Null terrain_path: - type: + type: - "string" - "null" description: "A specific path for the terrain file if using complex/interpolated terrain." default: Null bound_path: - type: + type: - "string" - "null" description: "A specific path to a MeshFunction file that stores the boundary IDs." @@ -187,42 +185,42 @@ properties: description: "Sets periodic boundary conditions along the y-direction (Not yet implemented)." default: False x_range: - type: + type: - "array" - "null" description: "The range of the domain in the streamwise direction, e.g., `[x_min, x_max]`." default: Null units: "meter" y_range: - type: + type: - "array" - "null" description: "The range of the domain in the spanwise direction, e.g., `[y_min, y_max]`." default: Null units: "meter" z_range: - type: + type: - "array" - "null" description: "The range of the domain in the vertical direction, e.g., `[z_min, z_max]`." default: Null units: "meter" nx: - type: + type: - "integer" - "null" description: "The integer number of nodes in the streamwise (x) direction." default: Null units: "dimensionless" ny: - type: + type: - "integer" - "null" description: "The integer number of nodes in the spanwise (y) direction." default: Null units: "dimensionless" nz: - type: + type: - "integer" - "null" description: "The integer number of nodes in the vertical (z) direction." @@ -237,34 +235,28 @@ properties: - "elliptic" - "squircular" - "stretch" + - "gmsh" center: - type: + type: - "array" - "null" description: "Center of the cylinder/circle domains, e.g., `[x_center, y_center]`." default: Null units: "meter" radius: - type: + type: - "number" - "null" description: "Radius for the cylinder/circle domains." default: Null units: "meter" nt: - type: + type: - "integer" - "null" description: "Integer number of nodes in the theta direction for cylinder/circle domains." default: Null units: "dimensionless" - res: - type: - - "number" - - "null" - description: "The characteristic cell length for cylinder/circle domains generated with the mshr mesh_type." - default: Null - units: "meter" interpolated: type: "boolean" description: "Lets you define a terrain_path to have complex domain." @@ -282,7 +274,7 @@ properties: - "array" - "null" description: "Center of the hill." - default: + default: - 0.0 - 0.0 units: "meter" @@ -292,21 +284,21 @@ properties: default: 0.0 units: "radians" amp: - type: + type: - "number" - "null" description: "Height of the hill." default: Null units: "meter" sigma_x: - type: + type: - "number" - "null" description: "Skew in x." default: Null units: "meter" sigma_y: - type: + type: - "number" - "null" description: "Skew in y." @@ -319,20 +311,20 @@ properties: intercept: type: "array" description: "Define a plane with z = mx(x-x0)+my(y-y0)+z0, where intercept = [x0, y0, z0]" - default: + default: - 0.0 - 0.0 - 0.0 units: "meter" mx: - type: + type: - "number" - "null" description: "x slope." default: Null units: "meter/meter" my: - type: + type: - "number" - "null" description: "y slope." @@ -347,7 +339,7 @@ properties: properties: type: default: Null - type: + type: - "string" - "null" enum: @@ -359,7 +351,7 @@ properties: description: "Type of wind farm, options are ``grid``, ``random``, ``imported``, or ``empty``. See :ref:`imported_windfarm`" path: default: Null - type: + type: - "string" - "null" description: "Location of imported wind farm." @@ -369,42 +361,42 @@ properties: description: "If true, then use matplotlib and show() the wind farm/chord profiles mid run." ex_x: default: Null - type: + type: - "array" - "null" description: "Extents of the farm in the x direction, e.g., [x_min, x_max]." units: "meter" ex_y: default: Null - type: + type: - "array" - "null" description: "Extents of the farm in the y direction, e.g., [y_min, y_max]." units: "meter" x_spacing: default: Null - type: + type: - "number" - "null" description: "x spacing between turbines." units: "meter" y_spacing: default: Null - type: + type: - "number" - "null" description: "y spacing between turbines." units: "meter" x_shear: default: Null - type: + type: - "number" - "null" description: "x offset between rows." units: "meter" y_shear: default: Null - type: + type: - "number" - "null" description: "y offset between columns." @@ -416,35 +408,35 @@ properties: units: "meter" grid_rows: default: Null - type: + type: - "integer" - "null" description: "Integer number of turbines in the y direction." units: "dimensionless" grid_cols: default: Null - type: + type: - "integer" - "null" description: "Integer number of turbines in the x direction." units: "dimensionless" jitter: default: 0.0 - type: + type: - "number" - "boolean" description: "Magnitude of random noise added to a gridded wind farm." units: "meter" numturbs: default: Null - type: + type: - "integer" - "null" description: "Total integer number of turbines." units: "dimensionless" seed: default: Null - type: + type: - "integer" - "null" description: "Seed to control/prescribe the randomness between runs." @@ -470,7 +462,7 @@ properties: properties: type: default: Null - type: + type: - "string" - "null" enum: @@ -486,35 +478,35 @@ properties: description: "Type of representation, options are ``2D_disk``, ``numpy_disk``, ``disk``, ``hybrid_disk``, ``line``, ``expr_disk``, ``dolfin_line``, or ``disabled``. See :ref:`turbine_representation`" HH: default: Null - type: + type: - "number" - "null" description: "Hub height." units: "meter" RD: default: Null - type: + type: - "number" - "null" description: "Rotor diameter." units: "meter" thickness: default: Null - type: + type: - "number" - "null" description: "Thickness of the actuator disk (usually 10% of RD)." units: "meter" yaw: default: Null - type: + type: - "number" - "null" description: "Yaw of the turbine relative to inflow angle." units: "radians" axial: default: Null - type: + type: - "number" - "null" description: "Axial induction value for actuator disks." @@ -529,14 +521,14 @@ properties: description: "Distribution of force along the radial direction of an actuator disk, options are ``constant``, ``sine``, or ``chord``." rpm: default: Null - type: + type: - "number" - "null" description: "Rotation specified for the alm method." units: "rotations/minute" read_turb_data: default: Null - type: + type: - "string" - "null" description: "Path to alm data." @@ -587,19 +579,19 @@ properties: units: "dimensionless" chord_override: default: Null - type: + type: - "string" - "null" description: "The path to a specific chord to use in CSV format, e.g., input_data/chord_base.csv." motion_file: default: Null - type: + type: - "string" - "null" description: "Location to the platform motion data." motion_type: default: Null - type: + type: - "array" - "string" - "null" @@ -621,7 +613,7 @@ properties: properties: warp_type: default: Null - type: + type: - "string" - "null" enum: @@ -631,21 +623,21 @@ properties: description: "Warping will shift the nodes along the z direction concentrating them near the ground, options are ``smooth`` or ``split``." warp_strength: default: Null - type: + type: - "number" - "null" description: "For smooth warps, how aggressively they are moved to the ground." units: "dimensionless" warp_percent: default: Null - type: + type: - "number" - "null" description: "For split warps, percentage moved below the warp_heigth." units: "dimensionless" warp_height: default: Null - type: + type: - "number" - "null" description: "For split warps, where the split happens." @@ -690,7 +682,7 @@ properties: units: "dimensionless" refine_custom: default: Null - type: + type: - "object" - "null" description: "Allows for a dictionary of custom refine commands. See :ref:`refine_custom`" @@ -707,7 +699,7 @@ properties: properties: type: default: Null - type: + type: - "string" - "null" description: "Type of function space, options are ``linear`` or ``taylor_hood``." @@ -738,7 +730,7 @@ properties: properties: vel_profile: default: Null - type: + type: - "string" - "null" description: "Inflow velocity profile, options are ``uniform``, ``power``, ``log``, or ``turbsim``." @@ -772,7 +764,7 @@ properties: units: "dimensionless" turbsim_path: default: Null - type: + type: - "string" - "null" description: "Location for the turbsim inflow data." @@ -789,49 +781,49 @@ properties: properties: east: default: Null - type: + type: - "string" - "null" description: "Positive x." north: default: Null - type: + type: - "string" - "null" description: "Positive y." west: default: Null - type: + type: - "string" - "null" description: "Negative x." south: default: Null - type: + type: - "string" - "null" description: "Negative y." bottom: default: Null - type: + type: - "string" - "null" description: "Negative z." top: default: Null - type: + type: - "string" - "null" description: "Positive z." inflow: default: Null - type: + type: - "string" - "null" description: "Inflow, used in cylinder/circle." outflow: default: Null - type: + type: - "string" - "null" description: "Outflow, used in cylinder/circle." @@ -841,25 +833,25 @@ properties: properties: inflow: default: Null - type: + type: - "array" - "null" description: "List of inflow bc names." no_slip: default: Null - type: + type: - "array" - "null" description: "List of no_slip bc names." free_slip: default: Null - type: + type: - "array" - "null" description: "List of free_slip bc names." no_stress: default: Null - type: + type: - "array" - "null" description: "List of no_stress bc names." @@ -873,7 +865,7 @@ properties: properties: type: default: Null - type: + type: - "string" - "null" description: "Type of model, options are ``stabilized``, ``steady``, ``taylor_hood``, ``iterative_steady``, or ``unsteady``." @@ -969,7 +961,7 @@ properties: description: "Include the final wind angle in the sweep." velocity_path: default: Null - type: + type: - "string" - "null" description: "Location of inflow velocities for multivelocity." @@ -1016,7 +1008,7 @@ properties: description: "Minimize or maximize." control_types: default: Null - type: + type: - "array" - "null" description: "Controls to optimize, list of: ``layout``, ``yaw``, ``axial``, ``chord``, ``lift``, or ``drag``." @@ -1051,7 +1043,7 @@ properties: units: "dimensionless" kwargs: default: Null - type: + type: - "string" - "null" description: "If constraint is based on an objective function with kwargs, set them here." @@ -1143,7 +1135,7 @@ properties: - write_floris_input properties: write_floris_input: - type: + type: - "array" - "null" default: Null