Add block force approximationts

environment.yaml

@@ -3,25 +3,24 @@ 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
   - pytest
   - 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 .

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

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"]
+  volume: [1e-9, "relative"]
+BoundaryManager:
+  avg_x:  [1e-3, "relative"]

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)
 

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=""):