Userdensity

windse/ProblemManager.py

@@ -594,9 +594,10 @@ def ComputeFunctional(self,inflow_angle):
         # Define fluid properties
         # FIXME: These should probably be set in params.yaml input filt
         # nu = 1/10000
-        rho = 1
+        rho = self.params["problem"]["density"] # User-defined change
         nu_c = Constant(self.viscosity, name="viscosity")
         rho_c = Constant(rho, name="rho")
+        MISSING_RHO_C = rho_c
 
         # Define time step size (this value is used only for step 1 if adaptive timestepping is used)
         # FIXME: change variable name to avoid confusion within dolfin adjoint
@@ -623,7 +624,7 @@ def ComputeFunctional(self,inflow_angle):
         self.u_k1 = Function(self.fs.V, name="u_k1")
         self.u_k2 = Function(self.fs.V, name="u_k2")
 
-        # Seed previous velocity fields with the chosen initial condition
+        # Seed previous velocity fields with the chosen initial condition #DEBUGING: Following lines are commented so that velocity isn't assigned to all nodes
         self.u_k.assign(self.bd.bc_velocity)
         self.u_k1.assign(self.bd.bc_velocity)
         self.u_k2.assign(self.bd.bc_velocity)
@@ -694,12 +695,19 @@ def ComputeFunctional(self,inflow_angle):
         #    + (nu_c+self.nu_T)*inner(grad(U_CN), grad(v))*dx \
         #    + dot(nabla_grad(self.p_k1), v)*dx \
         #    - dot(self.tf, v)*dx
+        # Update: density included with pressure term of NS
+        #F1 = (1.0/self.dt_c)*inner(u - self.u_k1, v)*dx \
+        #   + inner(dot(U_AB, nabla_grad(U_CN)), v)*dx \
+        #   + (nu_c+self.nu_T)*inner(grad(U_CN), grad(v))*dx \
+        #   + inner(grad(self.p_k1), v)*dx \
+        #   - dot(self.tf, v)*dx
 
         F1 = (1.0/self.dt_c)*inner(u - self.u_k1, v)*dx \
            + inner(dot(U_AB, nabla_grad(U_CN)), v)*dx \
            + (nu_c+self.nu_T)*inner(grad(U_CN), grad(v))*dx \
-           + inner(grad(self.p_k1), v)*dx \
-           - dot(self.tf, v)*dx
+           + (1.0/MISSING_RHO_C)*inner(grad(self.p_k1), v)*dx \
+           - (1.0/MISSING_RHO_C)*dot(self.tf, v)*dx
+
 
         self.a1 = lhs(F1)
         self.L1 = rhs(F1)
@@ -708,7 +716,8 @@ def ComputeFunctional(self,inflow_angle):
         # self.a2 = dot(nabla_grad(p), nabla_grad(q))*dx
         # self.L2 = dot(nabla_grad(self.p_k1), nabla_grad(q))*dx - (1.0/self.dt_c)*div(self.u_k)*q*dx
         self.a2 = inner(grad(p), grad(q))*dx
-        self.L2 = inner(grad(self.p_k1), grad(q))*dx - (1.0/self.dt_c)*div(self.u_k)*q*dx
+        #self.L2 = inner(grad(self.p_k1), grad(q))*dx - (1/self.dt_c)*div(self.u_k)*q*dx
+        self.L2 = inner(grad(self.p_k1), grad(q))*dx - (MISSING_RHO_C/self.dt_c)*div(self.u_k)*q*dx
 
         # phi = p - self.p_k
         # F2 = inner(grad(q), grad(phi))*dx - (1.0/self.dt_c)*div(u_k)*q*dx
@@ -719,7 +728,8 @@ def ComputeFunctional(self,inflow_angle):
         # self.a3 = dot(u, v)*dx
         # self.L3 = dot(self.u_k, v)*dx - self.dt_c*dot(nabla_grad(self.p_k - self.p_k1), v)*dx
         self.a3 = inner(u, v)*dx
-        self.L3 = inner(self.u_k, v)*dx - self.dt_c*inner(grad(self.p_k - self.p_k1), v)*dx
+        #self.L3 = inner(self.u_k, v)*dx - self.dt_c*inner(grad(self.p_k - self.p_k1), v)*dx
+        self.L3 = inner(self.u_k, v)*dx - (self.dt_c/MISSING_RHO_C)*inner(grad(self.p_k - self.p_k1), v)*dx
 
         # F3 = inner(u, v)*dx - inner(self.u_k, v)*dx + self.dt_c*inner(phi, v)*dx
         # self.a3 = lhs(F3)
@@ -744,4 +754,3 @@ def ComputeFunctional(self,inflow_angle):
 #             self.mcd[k] = Constant(cd[k])
 
 # # ================================================================
-

windse/SolverManager.py

@@ -1467,13 +1467,13 @@ def rot_z(theta):
         #================================================================
 
         # Set the density
-        rho = 1.0
+        rho = self.params["problem"]["density"]
 
         # Set the hub height
         hub_height = self.problem.farm.HH[0] # For a SWIFT turbine
 
         # Get the hub-height velocity
-        u_inf = 8.0
+        u_inf = 8.0 # DEBUGGING: This variable is set to a constant number. Is this right?
 
         # Set the rotational speed of the turbine
         RPM = 15.0
@@ -1482,7 +1482,7 @@ def rot_z(theta):
         yaw = self.problem.farm.yaw[0]
 
         # Set the number of blades in the turbine
-        num_blades = 3
+        num_blades = self.params["turbines"]["num_blades"]
 
         # Blade length (turbine radius)
         L = self.problem.farm.radius[0] # For a SWIFT turbine 27 m in diameter
@@ -1502,7 +1502,7 @@ def rot_z(theta):
         # print(self.problem.farm.radius)
 
         # Discretize each blade into separate nodes
-        num_actuator_nodes = 10
+        num_actuator_nodes = 10 # self.params["turbines"]["blade_segments"] # DEBUGGING: This variable is updated to be a user-defined parameter
 
         #================================================================
         # Set Derived Constants
@@ -1666,7 +1666,7 @@ def UpdateActuatorLineForceOld(self, simTime):
         coords = np.copy(coords[0::self.problem.dom.dim, :])
 
         # Set up the turbine geometry
-        num_blades = 3
+        num_blades = self.params["turbines"]["num_blades"]
         theta_vec = np.linspace(0, 2.0*np.pi, num_blades+1)
         theta_vec = theta_vec[0:num_blades]
 
@@ -1676,7 +1676,7 @@ def UpdateActuatorLineForceOld(self, simTime):
         cl = 1.0
         cd = 2.0
 
-        rho = 1
+        rho = self.params["problem"]["density"]
 
         u_inf = 8
 
@@ -1687,7 +1687,7 @@ def UpdateActuatorLineForceOld(self, simTime):
         c = L/20
 
         # Number of blade evaluation sections
-        num_actuator_nodes = 50
+        num_actuator_nodes = 50 # self.params["turbines"]["blade_segments"] # DEBUGGING: Hard-coded parameter changed to user-defined parameter
         rdim = np.linspace(0, L, num_actuator_nodes)
         zdim = 0.0 + np.zeros(num_actuator_nodes)
         xblade = np.vstack((np.zeros(num_actuator_nodes), rdim, zdim))

windse/input_schema.yaml

@@ -563,6 +563,11 @@ properties:
         type: "number"
         description: "The radius of the hub. If non-zero, actuator nodes will still be placed in the full range [0, rotor_radius], but the lift/drag properties in the range [0, hub_rad] will be modified to reflect the blade root."
         units: "meter"
+      num_blades:
+        default: 3
+        type: "integer"
+        description: "Integer number of blades for each turbine"
+        units: "dimensionless"
       chord_perturb:
         default: 0.0
         type: "number"
@@ -881,6 +886,11 @@ properties:
         type: "number"
         description: "Kinematic viscosity."
         units: "meter^2/second"
+      density:
+        default: 1.0
+        type: "number"
+        description: "Fluid density"
+        units: "kg/meter^3"
       lmax:
         default: 15
         type: "number"

windse/turbine_types/ActuatorDisk.py

@@ -202,7 +202,7 @@ def _setup_chord_force(self):
         if self.mchord is not None:
             if self.chord is not None:
                 if self.blade_segments == "computed":
-                    self.num_actuator_nodes = 10 ##### FIX THIS ####
+                    self.num_actuator_nodes = 10 ##### FIX THIS #### # DEBUGGING: self.num_actuator_nodes = int(2.0*self.radius/self.gaussian_width) # FromLine 150 of ActuatorLine.py
                     self.blade_segments = self.num_actuator_nodes
                 else:
                     self.num_actuator_nodes = self.blade_segments

windse/turbine_types/ActuatorLine.py

@@ -133,7 +133,7 @@ def init_constant_alm_terms(self, fs):
         #================================================================
 
         # Set the number of blades in the turbine
-        self.num_blades = 3
+        self.num_blades = self.params["turbines"]["num_blades"]
 
         # Set the spacing pf each blade
         self.theta_0_vec = np.linspace(0.0, 2.0*np.pi, self.num_blades, endpoint = False)
@@ -548,7 +548,7 @@ def get_angle_between_vectors(a, b, n):
                     tip_loss_fac[k] = 1.0
 
                 else:
-                    loss_exponent = 3.0/2.0*(self.radius-self.rdim[k])/(self.rdim[k]*np.sin(aoa))
+                    loss_exponent = self.num_blades/2.0*(self.radius-self.rdim[k])/(self.rdim[k]*np.sin(aoa)) #DEBUGGING: 3.0 is changed to self.num_blades for adjusted tip-loss factor
                     acos_arg = np.exp(-loss_exponent)
                     acos_arg = np.clip(acos_arg, -1.0, 1.0)
                     tip_loss_fac[k] = 2.0/np.pi*np.arccos(acos_arg)
@@ -676,11 +676,14 @@ def build_actuator_lines(self, u, inflow_angle, fs, dfd=None):
             else:
                 if self.lookup_cl is not None:
                     cl, cd, tip_loss_fac, aoa = self.lookup_lift_and_drag(u_rel, twist, self.blade_unit_vec[blade_id])
-
+            # self.fprint("For blade_id = {4}, the following values come from lookup_lift_and_drag: cl = {0}, cd = {1}, tip_loss_fac = {2}, aoa = {3} \n".format(cl,cd,tip_loss_fac,aoa,blade_id))
             # Calculate the lift and drag forces using the relative velocity magnitude
-            rho = 1.0
+            rho = self.params["problem"]["density"]
             lift = tip_loss_fac*(0.5*cl*rho*chord*self.w*u_rel_mag**2)
             drag = tip_loss_fac*(0.5*cd*rho*chord*self.w*u_rel_mag**2)
+
+            #self.fprint("For blade_id = {0} with rho = {1}, lift = {2}".format(blade_id,rho,lift)) #DEBUGGING
+            #self.fprint("drag = {}".format(drag)) #DEBUGGING
 
             # Tile the blade coordinates for every mesh point, [numGridPts*ndim x problem.num_actuator_nodes]
             blade_pos_full = np.tile(self.blade_pos[blade_id], (np.shape(self.coords)[0], 1))
@@ -726,19 +729,23 @@ def build_actuator_lines(self, u, inflow_angle, fs, dfd=None):
                 elif dfd == 'chord':
                     for j in range(self.ndim):
                         dfd_chord[j::self.ndim, k] += vector_nodal_lift[:, j] + vector_nodal_drag[:, j]
-
+                #self.fprint("Blade_id = {} \n".format(blade_id)) # DEBUGGING
                 # Compute the total force vector [x, y, z] at a single actuator node
                 actuator_lift = lift[k]*lift_unit_vec
+                #self.fprint("actuator_lift = {} \n".format(actuator_lift)) #DEBUGGING
                 actuator_drag = drag[k]*drag_unit_vec
+                #self.fprint("actuator_drag = {} \n".format(actuator_drag)) #DEBUGGING
 
                 # Note: since this will be used to define the force (torque) from fluid -> blade
                 # we reverse the direction that otherwise gives the turbine force from blade -> fluid
                 actuator_force = -(actuator_lift + actuator_drag)
+                #self.fprint("actuator_force = {} \n".format(actuator_force)) #DEBUGGING
                 # actuator_force = -(actuator_lift - actuator_drag)
 
                 # Find the component in the direction tangential to the blade
                 tangential_actuator_force = np.dot(actuator_force, self.blade_unit_vec[blade_id][:, 2])
 
+                #self.fprint("tangential_actuator_force = {} \n".format(tangential_actuator_force)) #DEBUGGING
                 rotor_plane_force = np.dot(actuator_force, self.blade_unit_vec[blade_id])
                 # fx.write('%.5f, ' % (rotor_plane_force[0]))
                 # fy.write('%.5f, ' % (rotor_plane_force[1]))
@@ -767,6 +774,7 @@ def build_actuator_lines(self, u, inflow_angle, fs, dfd=None):
             data = np.array(along_blade_quantities[save_val])
 
             nn = self.num_blades*self.num_actuator_nodes
+            # print(data) # DEBUGGING
 
             if np.size(data) == nn:
                 data = data.reshape(1, nn)
@@ -776,6 +784,7 @@ def build_actuator_lines(self, u, inflow_angle, fs, dfd=None):
 
             global_save_data = np.zeros((self.params.num_procs, nn))
             self.params.comm.Allgather(data, global_save_data)
+
             data = np.nanmean(data, axis=0)
 
             if self.first_call_to_alm:
@@ -864,17 +873,24 @@ def build_actuator_lines(self, u, inflow_angle, fs, dfd=None):
 
     def finalize_rotor_torque(self):
 
+        # self.fprint("num_procs = {} \n".format(self.params.num_procs)) #DEBUGGING
         data_in_torque = np.zeros(self.params.num_procs)
+        #self.fprint("data_in_torque = {} \n".format(data_in_torque)) #DEBUGGING
         self.params.comm.Gather(self.rotor_torque, data_in_torque, root=0)
+        #self.fprint("Gathered data_in_torque = {} \n".format(data_in_torque)) #DEBUGGING
 
         data_in_torque_count = np.zeros(self.params.num_procs, dtype=int)
+        #self.fprint("data_in_torque_count = {} \n".format(data_in_torque_count)) #DEBUGGING
         self.params.comm.Gather(self.rotor_torque_count, data_in_torque_count, root=0)
+        #self.fprint("Gathered data_in_torque_count = {} \n".format(data_in_torque_count)) #DEBUGGING
 
         if self.params.rank == 0:
             rotor_torque_sum = np.sum(data_in_torque)
+            #self.fprint("rotor_torque_sum = {} \n".format(rotor_torque_sum)) #DEBUG
             rotor_torque_count_sum = np.sum(data_in_torque_count)
-
-            # This removes the possibility of a power being doubled or tripled
+            #self.fprint("rotor_torque_count_sum = {} \n".format(rotor_torque_count_sum)) #DEBUG
+
+           # This removes the possibility of a power being doubled or tripled
             # if multiple ranks include this turbine and therefore calculate a torque
             self.rotor_torque = rotor_torque_sum/rotor_torque_count_sum
 
@@ -913,6 +929,7 @@ def turbine_force(self, u, inflow_angle, fs, **kwargs):
 
         # Call the function to build the complete mpi_u_fluid array
         self.get_u_fluid_at_alm_nodes(u)
+        #self.fprint("mpi_u_fluid = {}".format(self.mpi_u_fluid)) # DEBUGGING
 
         # Call the ALM function for this turbine
         if update:

windse/turbine_types/ActuatorLineDolfin.py

@@ -79,6 +79,8 @@ def load_parameters(self):
         self.motion_type = self.params["turbines"]["motion_type"]
         self.use_gauss_vel_probe = self.params["turbines"]["use_gauss_vel_probe"]
         self.use_ap_linear_interp = self.params["turbines"]["use_ap_linear_interp"]
+        self.num_blades = self.params["turbines"]["num_blades"]
+        self.density = self.params["problem"]["density"]
 
     def compute_parameters(self):
         pass
@@ -172,8 +174,8 @@ def calc_platform_motion(self,T):
 
     def init_alm_calc_terms(self):
         # compute turbine radius
-        self.rho = 1.0
-        self.num_blades = 3
+        self.rho = self.density
+        self.num_blades = self.num_blades
         self.radius = 0.5*self.RD
         self.angular_velocity = 2.0*pi*self.rpm/60.0
 
@@ -594,7 +596,7 @@ def calculate_tip_loss(self, rdim, aoa):
                 tip_loss_fac = 1.0
 
             else:
-                loss_exponent = 3.0/2.0*(self.radius-rdim)/(rdim*sin(aoa))
+                loss_exponent = self.num_blades/2.0*(self.radius-rdim)/(rdim*sin(aoa)) # DEBUGGING: 3.0 is changed to self.num_blades
                 acos_arg = exp(-loss_exponent)
         #         acos_arg = np.clip(acos_arg, -1.0, 1.0)
                 tip_loss_fac = 2.0/pi*acos(acos_arg)
@@ -647,6 +649,7 @@ def build_actuator_node(self, u_k, x_0, x_0_pre_motion, n_0, blade_id, actuator_
         w = self.w[actuator_id]
         twist = self.mtwist[actuator_id]
         chord = self.mchord[actuator_id]
+        self.rho = self.density # Density is redefined for calculation of lift and drag
 
         # Build a spherical Gaussian with width eps about point x_0
         self.gauss_kernel[blade_id][actuator_id] = self.build_gauss_kernel_at_point(x_0)