diff --git a/.github/workflows/test_pvade.yaml b/.github/workflows/test_pvade.yaml index dcb979b2..f9e3b437 100644 --- a/.github/workflows/test_pvade.yaml +++ b/.github/workflows/test_pvade.yaml @@ -33,6 +33,8 @@ jobs: shell: bash -l {0} run: PYTHONPATH=. pytest -sv pvade/tests/ - name: test all inputs + env: + OMPI_MCA_regx: "naive" shell: bash -l {0} run: pytest -sv test_all_inputs.py diff --git a/examples/cylinderflow.py b/examples/cylinderflow.py index 8aae4991..e1b9600f 100644 --- a/examples/cylinderflow.py +++ b/examples/cylinderflow.py @@ -45,7 +45,6 @@ from dolfinx.io import VTXWriter, gmshio, XDMFFile from dolfinx.mesh import locate_entities_boundary - #################################################### # # # MESH PARAMETERS # diff --git a/examples/poissoneq.py b/examples/poissoneq.py index f2968894..96b55728 100644 --- a/examples/poissoneq.py +++ b/examples/poissoneq.py @@ -10,7 +10,6 @@ from timeit import default_timer as timer from dolfinx.common import TimingType, list_timings - start = timer() elems = int(sys.argv[1]) # nelems # Create mesh and define function space diff --git a/examples/synthetic_turbulent_inflow/generate_turbulent_inflow_h5_file.ipynb b/examples/synthetic_turbulent_inflow/generate_turbulent_inflow_h5_file.ipynb index be392721..1a22d114 100644 --- a/examples/synthetic_turbulent_inflow/generate_turbulent_inflow_h5_file.ipynb +++ b/examples/synthetic_turbulent_inflow/generate_turbulent_inflow_h5_file.ipynb @@ -24,12 +24,18 @@ "import pandas as pd # need this to load our data from the csv files\n", "\n", "import sys\n", + "\n", "# NOTE: to run this, you will need to replace the path below #\n", - "sys.path.append('/Users/bstanisl/repos/pyconturb/pyconturb') # path to the pyconturb repository location on your local machine\n", + "sys.path.append(\n", + " \"/Users/bstanisl/repos/pyconturb/pyconturb\"\n", + ") # path to the pyconturb repository location on your local machine\n", "from pyconturb import gen_turb, gen_spat_grid # generate turbulence, useful helper\n", "from pyconturb.sig_models import iec_sig # IEC 61400-1 turbulence std dev\n", "from pyconturb.spectral_models import kaimal_spectrum # Kaimal spectrum\n", - "from pyconturb.wind_profiles import constant_profile, power_profile # wind-speed profile functions\n", + "from pyconturb.wind_profiles import (\n", + " constant_profile,\n", + " power_profile,\n", + ") # wind-speed profile functions\n", "\n", "from _nb_utils import plot_slice\n", "import h5py" @@ -217,11 +223,13 @@ "source": [ "# generate spatial gridpoints dataframe\n", "ny = 80\n", - "nz = 80 #40\n", + "nz = 80 # 40\n", "y = np.linspace(-10.0, 10.0, ny)\n", - "z = np.linspace(0.00001, 20.0, nz) \n", + "z = np.linspace(0.00001, 20.0, nz)\n", "\n", - "spat_df = gen_spat_grid(y, z) # if `comps` not passed in, assumes all 3 components are wanted\n", + "spat_df = gen_spat_grid(\n", + " y, z\n", + ") # if `comps` not passed in, assumes all 3 components are wanted\n", "spat_df.head() # look at the first few rows" ] }, @@ -233,9 +241,9 @@ "outputs": [], "source": [ "# generate time vector\n", - "t_final = 1.0 #10.0 [s]\n", - "dt = 0.01 # [s]\n", - "t_steps = int(t_final/dt)\n", + "t_final = 1.0 # 10.0 [s]\n", + "dt = 0.01 # [s]\n", + "t_steps = int(t_final / dt)\n", "time = np.linspace(0, t_final, t_steps)" ] }, @@ -483,16 +491,30 @@ "source": [ "# reshape to 3D array and visualize first timestep\n", "data = {}\n", - "data['u'] = turb_df.filter(regex='u').values.reshape(len(turb_df),y.size,z.size).transpose((0, 2, 1))\n", - "data['v'] = turb_df.filter(regex='v').values.reshape(len(turb_df),y.size,z.size).transpose((0, 2, 1))\n", - "data['w'] = turb_df.filter(regex='w').values.reshape(len(turb_df),y.size,z.size).transpose((0, 2, 1))\n", + "data[\"u\"] = (\n", + " turb_df.filter(regex=\"u\")\n", + " .values.reshape(len(turb_df), y.size, z.size)\n", + " .transpose((0, 2, 1))\n", + ")\n", + "data[\"v\"] = (\n", + " turb_df.filter(regex=\"v\")\n", + " .values.reshape(len(turb_df), y.size, z.size)\n", + " .transpose((0, 2, 1))\n", + ")\n", + "data[\"w\"] = (\n", + " turb_df.filter(regex=\"w\")\n", + " .values.reshape(len(turb_df), y.size, z.size)\n", + " .transpose((0, 2, 1))\n", + ")\n", "\n", "fig, ax = plt.subplots()\n", - "plt.imshow(data['u'][0,:,:], # imshow requires nz-ny slice\n", - " origin='lower', # smallest y-z in lower left, not upper left\n", - " extent=[y[0], y[-1], z[0], z[-1]], # lateral and vertical limits\n", - " interpolation='bilinear',\n", - " cmap='coolwarm') # image smoothing\n", + "plt.imshow(\n", + " data[\"u\"][0, :, :], # imshow requires nz-ny slice\n", + " origin=\"lower\", # smallest y-z in lower left, not upper left\n", + " extent=[y[0], y[-1], z[0], z[-1]], # lateral and vertical limits\n", + " interpolation=\"bilinear\",\n", + " cmap=\"coolwarm\",\n", + ") # image smoothing\n", "plt.colorbar()" ] }, @@ -507,24 +529,24 @@ "h5_filename = \"pct_turb_ny{}_nz{}_unconstrained_{}s_dt{}_uref{}.h5\".format(\n", " ny, nz, t_final, dt, u_ref\n", ")\n", - "with h5py.File(\"../../input/\"+h5_filename, \"w\") as fp:\n", + "with h5py.File(\"../../input/\" + h5_filename, \"w\") as fp:\n", " fp.create_dataset(\"time_index\", shape=(t_steps,))\n", " fp[\"time_index\"][:] = time\n", - " \n", + "\n", " fp.create_dataset(\"y_coordinates\", shape=(ny,))\n", " fp[\"y_coordinates\"][:] = y\n", - " \n", + "\n", " fp.create_dataset(\"z_coordinates\", shape=(nz,))\n", " fp[\"z_coordinates\"][:] = z\n", - " \n", + "\n", " fp.create_dataset(\"u\", shape=(t_steps, nz, ny))\n", - " fp[\"u\"][:] = data['u'][:]\n", - " \n", + " fp[\"u\"][:] = data[\"u\"][:]\n", + "\n", " fp.create_dataset(\"v\", shape=(t_steps, nz, ny))\n", - " fp[\"v\"][:] = data['v'][:]\n", - " \n", + " fp[\"v\"][:] = data[\"v\"][:]\n", + "\n", " fp.create_dataset(\"w\", shape=(t_steps, nz, ny))\n", - " fp[\"w\"][:] = data['w'][:]" + " fp[\"w\"][:] = data[\"w\"][:]" ] }, { @@ -551,7 +573,7 @@ "source": [ "# read h5 file\n", "# Note: Replace the path below with the path on your local machine to the .h5 file that you want to read\n", - "fname = '/Users/bstanisl/Documents/repos/PVade/input/pct_turb_ny80_nz80_unconstrained_1.0s_dt0.01_uref20.h5'\n", + "fname = \"/Users/bstanisl/Documents/repos/PVade/input/pct_turb_ny80_nz80_unconstrained_1.0s_dt0.01_uref20.h5\"\n", "\n", "# Open the file (read-only)\n", "with h5py.File(fname, \"r\") as f:\n", diff --git a/input/duramat_case_study.yaml b/input/duramat_case_study.yaml index fb5939e9..646b0a1d 100644 --- a/input/duramat_case_study.yaml +++ b/input/duramat_case_study.yaml @@ -23,6 +23,13 @@ pv_array: elevation: 2.1 tracker_angle: 0.0 span_fixation_pts: [13.2] + # torque_tube_separation: 0.2 # gap between panel and tube center + # torque_tube_outer_radius: 0.1 # radius of the torque tube + # torque_tube_inner_radius: 0.09 # radius of the torque tube + modules_per_span: 1 + fixed_location: 1 # fixed location of the panel along the span (from 0 (fixed at left) to modules_per_span (fixed at right)) + block_chord_div_by_panel_chord: 0.02 + solver: dt: 0.01 t_final: 20.0 @@ -58,4 +65,7 @@ structure: bc_list: [] motor_connection: true tube_connection: true - beta_relaxation: 0.5 \ No newline at end of file + beta_relaxation: 0.5 + elasticity_modulus_tube: 2.0e+11 + poissons_ratio_tube: 0.3 + rho_tube: 7800.0 \ No newline at end of file diff --git a/input/flag2d.yaml b/input/flag2d.yaml index 561a698e..307d36db 100644 --- a/input/flag2d.yaml +++ b/input/flag2d.yaml @@ -52,7 +52,7 @@ structure: body_force_x: 0 body_force_y: 0 body_force_z: 0 #100 - bc_list: ["left"] + bc_list: ["panel_left"] motor_connection: False tube_connection: False beta_relaxation: 0.005 diff --git a/input/heated_panels2d.yaml b/input/heated_panels2d.yaml index 468634ff..dfc7c946 100644 --- a/input/heated_panels2d.yaml +++ b/input/heated_panels2d.yaml @@ -59,7 +59,7 @@ structure: body_force_x: 0 body_force_y: -1 body_force_z: 0 #100 - bc_list: ["left"] + bc_list: ["panel_left"] motor_connection: False tube_connection: False beta_relaxation: 0.005 diff --git a/input/inflow_input.yaml b/input/inflow_input.yaml index 539e87f0..0646adfa 100644 --- a/input/inflow_input.yaml +++ b/input/inflow_input.yaml @@ -56,7 +56,7 @@ structure: body_force_x: 0 body_force_y: 0 body_force_z: 0 #100 - bc_list: ["left"] + bc_list: ["panel_left"] motor_connection: False tube_connection: False beta_relaxation: 0.005 diff --git a/input/panels2d.yaml b/input/panels2d.yaml index 7482b212..57d80cd9 100644 --- a/input/panels2d.yaml +++ b/input/panels2d.yaml @@ -46,7 +46,7 @@ structure: body_force_x: 0 body_force_y: -1 body_force_z: 0 #100 - bc_list: ["top"] + bc_list: [panel_top] motor_connection: False tube_connection: False beta_relaxation: 0.005 diff --git a/input/panels3d.yaml b/input/panels3d.yaml index a66e1c53..a6159f10 100644 --- a/input/panels3d.yaml +++ b/input/panels3d.yaml @@ -56,5 +56,5 @@ structure: body_force_x: 0 body_force_y: 0 body_force_z: -1 #100 - bc_list: [left ] + bc_list: [panel_left] tube_connection: False diff --git a/pvade/IO/input_schema.yaml b/pvade/IO/input_schema.yaml index 4aefbfa4..b389509a 100644 --- a/pvade/IO/input_schema.yaml +++ b/pvade/IO/input_schema.yaml @@ -136,7 +136,7 @@ properties: minimum: 0.0 maximum: 10.0 type: "number" - description: "The vertical distance between the center of the panel and the ground." + description: "The vertical distance between the center of the torque tube and the ground." units: "meter" stream_spacing: default: 7.0 @@ -154,7 +154,7 @@ properties: units: "meter" span_fixation_pts: default: 3.5 - minimum: 1.0 + minimum: 0.2 # maximum: 32.0 type: - "number" @@ -198,6 +198,39 @@ properties: be unrolled such that tracking_angle_0 is applied to the panel at (x_min, y_min), tracking_angle_1 is applied to the panel at (x_min + x_spacing, y_min), etc., i.e., x-major direction. If argument is a single value, that tracking angle will be applied to the every panel." units: "degree" + torque_tube_separation: + default: 0.0 + minimum: 0.0 + maximum: 1.0 + type: "number" + description: "The distance between the bottom of the panel structure and the centerline of the torque tube, in meters." + torque_tube_outer_radius: + default: 0.0 + minimum: 0.0 + maximum: 1.0 + type: "number" + description: "The radius of the torque tube, in meters." + torque_tube_inner_radius: + default: 0.0 + minimum: 0.0 + maximum: 1.0 + type: "number" + description: "The radius of the torque tube, in meters." + modules_per_span: + default: 1 + minimum: 1 + type: "integer" + description: "The number of modules/panels per every panel_span distance. This can be thought of as the number of modules/panels per each table, counted in the spanwise direction only (i.e., 2-in-portait systems should report only the number in the spanwise direction, vs 2x that value)." + fixed_location: + default: 5 + type: "integer" + description: "If an integer between 0 and modules_per_span, indicates the fixed location" + block_chord_div_by_panel_chord: + default: 0.1 + minimum: 0.01 + type: "number" + description: "block length or width divided by panel chord" + solver: additionalProperties: false type: "object" @@ -617,7 +650,7 @@ properties: type: "boolean" description: "Controls whether or not a boundary condition is applied along the lines defined by the spanwise fixation points which run in the streamwise direction and divide the panel into multiple rectangles in the spanwise direction." bc_list: - default: [front,back] + default: []#[panel_front,panel_back] # type: "list" description: "location for clamped boundary conditions" dt: @@ -648,6 +681,27 @@ properties: type: "number" description: "poissons ratio of the panel structure." units: "None" + elasticity_modulus_tube: + default: 200.0e+9 + minimum: 1.0e+2 + maximum: 500.0e+9 + type: "number" + description: "The effective Young's modulus of the torque tube." + units: "Pascal" + poissons_ratio_tube: + default: 0.3 + minimum: 0.1 + maximum: 0.9 + type: "number" + description: "poissons ratio of the torque tube." + units: "None" + rho_tube: + default: 1.0 + # minimum: 0.001 + # maximum: 1000.0 + type: "number" + description: "The density of the structure." + units: "kg/meter^3" body_force_x: default: 100 type: "number" diff --git a/pvade/fluid/FlowManager.py b/pvade/fluid/FlowManager.py index bd2c5c41..640e9713 100644 --- a/pvade/fluid/FlowManager.py +++ b/pvade/fluid/FlowManager.py @@ -165,6 +165,19 @@ def build_forms(self, domain, params): """ + # initialize total torque on each panel to zero + for panel_id in range( + int(params.pv_array.stream_rows * params.pv_array.span_rows) + ): + attr_name = f"total_torque_panel_{panel_id:.0f}" + setattr(self, attr_name, dolfinx.fem.Constant(domain.fluid.msh, 0.0)) + + for module_id in range(params.pv_array.modules_per_span + 1): + attr_d_name = ( + f"double_integral_total_torque_panel_{panel_id:.0f}_{module_id:.0f}" + ) + setattr(self, attr_d_name, dolfinx.fem.Constant(domain.fluid.msh, 0.0)) + # Define fluid properties if self.ndim == 2: self.dpdx = dolfinx.fem.Constant(domain.fluid.msh, (0.0, 0.0)) @@ -242,6 +255,50 @@ def build_forms(self, domain, params): self.p_k = dolfinx.fem.Function(self.Q, name="pressure") self.p_k1 = dolfinx.fem.Function(self.Q) + # define function to save x, y, z coordinates + self.spatial_X_coords = dolfinx.fem.Function(self.Q, name="X_coords") + self.spatial_Y_coords = dolfinx.fem.Function(self.Q, name="Y_coords") + self.spatial_Z_coords = dolfinx.fem.Function(self.Q, name="Z_coords") + + class FillFunctionWithXCoords: + def __init__(self): + pass + + def __call__(self, x): + coords = np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType) + + coords[0] = x[0] + + return coords + + self.spatial_X_coords.interpolate(FillFunctionWithXCoords()) + + class FillFunctionWithYCoords: + def __init__(self): + pass + + def __call__(self, x): + coords = np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType) + + coords[0] = x[1] + + return coords + + self.spatial_Y_coords.interpolate(FillFunctionWithYCoords()) + + class FillFunctionWithZCoords: + def __init__(self): + pass + + def __call__(self, x): + coords = np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType) + + coords[0] = x[2] + + return coords + + self.spatial_Z_coords.interpolate(FillFunctionWithZCoords()) + # initial conditions if params.fluid.initialize_with_inflow_bc: # self.inflow_profile = dolfinx.fem.Function(self.V) @@ -535,6 +592,37 @@ def _all_interior_surfaces(x): "ds", domain=domain.fluid.msh, subdomain_data=domain.fluid.facet_tags ) + # def compute_total_torques_on_each_panel(self, domain, params): + + # ds_fluid = ufl.Measure( + # "ds", domain=domain.fluid.msh, subdomain_data=domain.fluid.facet_tags + # ) + + # for panel_id in range( + # int(params.pv_array.stream_rows * params.pv_array.span_rows) + # ): + # total_torque = 0 + + # for module_id in range(params.pv_array.modules_per_span): + + # total_torque += dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(self.spatial_X_coords, self.traction[2]) * ds_fluid( + # domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"])) + # ) + # total_torque += dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(self.spatial_X_coords, self.traction[2]) * ds_fluid( + # domain.domain_markers[f"panel_top_{panel_id:.0f}_{module_id:.0f}"]["idx"])) + # ) + # total_torque += dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(self.spatial_Z_coords, self.traction[0]) * ds_fluid( + # domain.domain_markers[f"panel_top_{panel_id:.0f}_{module_id:.0f}"]["idx"])) + # ) + # total_torque += dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(self.spatial_Z_coords, self.traction[0]) * ds_fluid( + # domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"])) + # ) + + # attr_name = f"total_torque_panel_{panel_id:.0f}" + + # setattr(self, attr_name, total_torque) + + # define integrated force forms self.integrated_force_x_form = [] self.integrated_force_y_form = [] self.integrated_force_z_form = [] @@ -548,80 +636,446 @@ def _all_interior_surfaces(x): # print(f"loc = {loc}, idx = {idx}, s = {s}") self.integrated_force_x_form.append(0) + self.integrated_force_y_form.append(0) + self.integrated_force_z_form.append(0) + self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid( - domain.domain_markers[f"left_{panel_id:.0f}"]["idx"] + domain.domain_markers[f"panel_left_{panel_id:.0f}"]["idx"] ) self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid( - domain.domain_markers[f"top_{panel_id:.0f}"]["idx"] + domain.domain_markers[f"panel_right_{panel_id:.0f}"]["idx"] ) - self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid( - domain.domain_markers[f"right_{panel_id:.0f}"]["idx"] + self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid( + domain.domain_markers[f"panel_left_{panel_id:.0f}"]["idx"] ) - self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid( - domain.domain_markers[f"bottom_{panel_id:.0f}"]["idx"] + self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid( + domain.domain_markers[f"panel_right_{panel_id:.0f}"]["idx"] ) + if self.ndim == 3: self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid( - domain.domain_markers[f"front_{panel_id:.0f}"]["idx"] + domain.domain_markers[f"panel_front_{panel_id:.0f}"]["idx"] + ) + self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid( + domain.domain_markers[f"panel_back_{panel_id:.0f}"]["idx"] + ) + self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid( + domain.domain_markers[f"panel_front_{panel_id:.0f}"]["idx"] + ) + self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid( + domain.domain_markers[f"panel_back_{panel_id:.0f}"]["idx"] ) + self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid( + domain.domain_markers[f"panel_front_{panel_id:.0f}"]["idx"] + ) + self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid( + domain.domain_markers[f"panel_back_{panel_id:.0f}"]["idx"] + ) + self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid( + domain.domain_markers[f"panel_left_{panel_id:.0f}"]["idx"] + ) + self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid( + domain.domain_markers[f"panel_right_{panel_id:.0f}"]["idx"] + ) + for module_id in range(params.pv_array.modules_per_span): + self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid( - domain.domain_markers[f"back_{panel_id:.0f}"]["idx"] + domain.domain_markers[ + f"panel_bottom_{panel_id:.0f}_{module_id:.0f}" + ]["idx"] + ) + self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid( + domain.domain_markers[f"panel_top_{panel_id:.0f}_{module_id:.0f}"][ + "idx" + ] + ) + self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid( + domain.domain_markers[ + f"panel_bottom_{panel_id:.0f}_{module_id:.0f}" + ]["idx"] + ) + self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid( + domain.domain_markers[f"panel_top_{panel_id:.0f}_{module_id:.0f}"][ + "idx" + ] ) + if self.ndim == 3: + + self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid( + domain.domain_markers[ + f"panel_bottom_{panel_id:.0f}_{module_id:.0f}" + ]["idx"] + ) + self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid( + domain.domain_markers[ + f"panel_top_{panel_id:.0f}_{module_id:.0f}" + ]["idx"] + ) + self.integrated_force_x_form[-1] = dolfinx.fem.form( self.integrated_force_x_form[-1] ) + self.integrated_force_y_form[-1] = dolfinx.fem.form( + self.integrated_force_y_form[-1] + ) + self.integrated_force_z_form[-1] = dolfinx.fem.form( + self.integrated_force_z_form[-1] + ) - self.integrated_force_y_form.append(0) - self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid( - domain.domain_markers[f"left_{panel_id:.0f}"]["idx"] + def compute_double_integral_panel_torques(self, domain, params): + + ds_fluid = ufl.Measure( + "ds", domain=domain.fluid.msh, subdomain_data=domain.fluid.facet_tags + ) + + for panel_id in range( + int(params.pv_array.stream_rows * params.pv_array.span_rows) + ): + whole_top_surface_facet_index = ( + [] + ) # facet index of the whole panel top surface in a row + whole_bot_surface_facet_index = ( + [] + ) # facet index of the whole panel bottom surface in a row + + for module_id in range(params.pv_array.modules_per_span): + whole_top_surface_facet_index += domain.fluid.facet_tags.find( + domain.domain_markers[f"panel_top_{panel_id:.0f}_{module_id:.0f}"][ + "idx" + ] + ).tolist() + whole_bot_surface_facet_index += domain.fluid.facet_tags.find( + domain.domain_markers[ + f"panel_bottom_{panel_id:.0f}_{module_id:.0f}" + ]["idx"] + ).tolist() + + whole_top_surf_submesh, entity_map, vertex_map, geom_map = ( + dolfinx.mesh.create_submesh( + domain.fluid.msh, + 2, + np.array(whole_top_surface_facet_index, dtype=np.int32), + ) ) - self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid( - domain.domain_markers[f"top_{panel_id:.0f}"]["idx"] + whole_bot_surf_submesh, entity_map, vertex_map, geom_map = ( + dolfinx.mesh.create_submesh( + domain.fluid.msh, + 2, + np.array(whole_bot_surface_facet_index, dtype=np.int32), + ) ) - self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid( - domain.domain_markers[f"right_{panel_id:.0f}"]["idx"] + + spatial_X_coords_top = dolfinx.fem.Function( + dolfinx.fem.FunctionSpace(whole_top_surf_submesh, ("CG", 1)), + name="X_coords_top", ) - self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid( - domain.domain_markers[f"bottom_{panel_id:.0f}"]["idx"] + spatial_X_coords_bot = dolfinx.fem.Function( + dolfinx.fem.FunctionSpace(whole_bot_surf_submesh, ("CG", 1)), + name="X_coords_bot", ) - if self.ndim == 3: - self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid( - domain.domain_markers[f"front_{panel_id:.0f}"]["idx"] - ) - self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid( - domain.domain_markers[f"back_{panel_id:.0f}"]["idx"] - ) + spatial_Z_coords_top = dolfinx.fem.Function( + dolfinx.fem.FunctionSpace(whole_top_surf_submesh, ("CG", 1)), + name="Z_coords_top", + ) + spatial_Z_coords_bot = dolfinx.fem.Function( + dolfinx.fem.FunctionSpace(whole_bot_surf_submesh, ("CG", 1)), + name="Z_coords_bot", + ) + spatial_X_coords_top.interpolate(self.spatial_X_coords) + spatial_X_coords_bot.interpolate(self.spatial_X_coords) + spatial_Z_coords_top.interpolate(self.spatial_Z_coords) + spatial_Z_coords_bot.interpolate(self.spatial_Z_coords) + """ + the torque should be relative to central of tube, traction x times z (z=0 at center tube) + traction z times x, (x=0 at center tube) + for z, how we should put the tube center at zero? It is not minus elevation + """ - self.integrated_force_y_form[-1] = dolfinx.fem.form( - self.integrated_force_y_form[-1] + whole_top_submesh_function_space = dolfinx.fem.FunctionSpace( + whole_top_surf_submesh, ("CG", 1) + ) + whole_bot_submesh_function_space = dolfinx.fem.FunctionSpace( + whole_bot_surf_submesh, ("CG", 1) ) - self.integrated_force_z_form.append(0) - if self.ndim == 3: - self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid( - domain.domain_markers[f"left_{panel_id:.0f}"]["idx"] + traction_z_array_top = dolfinx.fem.Function( + whole_top_submesh_function_space, name="traction_z_top" + ) + traction_x_array_top = dolfinx.fem.Function( + whole_top_submesh_function_space, name="traction_x_top" + ) + traction_z_array_bot = dolfinx.fem.Function( + whole_bot_submesh_function_space, name="traction_z_bot" + ) + traction_x_array_bot = dolfinx.fem.Function( + whole_bot_submesh_function_space, name="traction_x_bot" + ) + + traction_z_fluid = dolfinx.fem.Function(self.Q, name="traction_z_fluid") + traction_x_fluid = dolfinx.fem.Function(self.Q, name="traction_x_fluid") + + # traction is expression on whole fluid mesh, we need to + # project it to a function defined on whole fluid mesh, then + # interpolate it to the top and bot submesh function space. + + u_inter_fluid = ufl.TrialFunction(self.Q) + v_inter_fluid = ufl.TestFunction(self.Q) + + a_inter_fluid = ufl.inner(u_inter_fluid, v_inter_fluid) * ufl.ds + L_inter_fluid_z = ufl.inner(self.traction[2], v_inter_fluid) * ufl.ds + L_inter_fluid_x = ufl.inner(self.traction[0], v_inter_fluid) * ufl.ds + + # Assemble and solve + A_inter_fluid = dolfinx.fem.petsc.assemble_matrix( + dolfinx.fem.form(a_inter_fluid) + ) + A_inter_fluid.assemble() + + b_inter_fluid_z = dolfinx.fem.petsc.assemble_vector( + dolfinx.fem.form(L_inter_fluid_z) + ) + b_inter_fluid_x = dolfinx.fem.petsc.assemble_vector( + dolfinx.fem.form(L_inter_fluid_x) + ) + + solver = PETSc.KSP().create(traction_z_array_top.function_space.mesh.comm) + solver.setOperators(A_inter_fluid) + solver.setType("cg") + solver.getPC().setType("jacobi") + solver.setFromOptions() + solver.solve(b_inter_fluid_z, traction_z_fluid.vector) + traction_z_fluid.x.scatter_forward() + + solver = PETSc.KSP().create(traction_x_array_top.function_space.mesh.comm) + solver.setOperators(A_inter_fluid) + solver.setType("cg") + solver.getPC().setType("jacobi") + solver.setFromOptions() + solver.solve(b_inter_fluid_x, traction_x_fluid.vector) + traction_x_fluid.x.scatter_forward() + + traction_z_array_top.interpolate(traction_z_fluid) + traction_x_array_top.interpolate(traction_x_fluid) + traction_z_array_bot.interpolate(traction_z_fluid) + traction_x_array_bot.interpolate(traction_x_fluid) + + # traction_z_array_top.x.array[:] = 10*np.sin(params.pv_array.tracker_angle[0]*np.pi/180) + # traction_z_array_bot.x.array[:] = 10*np.sin(params.pv_array.tracker_angle[0]*np.pi/180) + # traction_x_array_top.x.array[:] = 10*np.cos(params.pv_array.tracker_angle[0]*np.pi/180) + # traction_x_array_bot.x.array[:] = 10*np.cos(params.pv_array.tracker_angle[0]*np.pi/180) + torque_function_on_panels_top = dolfinx.fem.Function( + whole_top_submesh_function_space + ) + torque_function_on_panels_top.x.array[:] = ( + spatial_X_coords_top.x.array[:] + - np.linspace( + 0.0, + params.pv_array.stream_spacing * (params.pv_array.stream_rows - 1), + params.pv_array.stream_rows, + )[panel_id] + ) * traction_z_array_top.x.array[:] - ( + spatial_Z_coords_top.x.array[:] - params.pv_array.elevation + ) * traction_x_array_top.x.array[ + : + ] + + torque_function_on_panels_bot = dolfinx.fem.Function( + whole_bot_submesh_function_space + ) + torque_function_on_panels_bot.x.array[:] = ( + spatial_X_coords_bot.x.array[:] + - np.linspace( + 0.0, + params.pv_array.stream_spacing * (params.pv_array.stream_rows - 1), + params.pv_array.stream_rows, + )[panel_id] + ) * traction_z_array_bot.x.array[:] - ( + spatial_Z_coords_bot.x.array[:] - params.pv_array.elevation + ) * traction_x_array_bot.x.array[ + : + ] + dx_top_whole = ufl.Measure("dx", domain=whole_top_surf_submesh) + dx_bot_whole = ufl.Measure("dx", domain=whole_bot_surf_submesh) + + total_torque_on_panel_array = dolfinx.fem.assemble_scalar( + dolfinx.fem.form(torque_function_on_panels_top * dx_top_whole) + ) + dolfinx.fem.assemble_scalar( + dolfinx.fem.form(torque_function_on_panels_bot * dx_bot_whole) + ) + + attr_name = f"total_torque_panel_{panel_id:.0f}" + total_torque_constant = getattr(self, attr_name) + + total_torque_constant.value = total_torque_on_panel_array + + coords_top = whole_top_submesh_function_space.tabulate_dof_coordinates() + coords_bot = whole_bot_submesh_function_space.tabulate_dof_coordinates() + + top_integrated_func_along_y = np.zeros_like( + torque_function_on_panels_top.x.array[:] + ) + bot_integrated_func_along_y = np.zeros_like( + torque_function_on_panels_bot.x.array[:] + ) + + # in PVade, the front has the minimum y coordinates, if we integral from large y to small y, it is negative, + # to keep it consistant as the theoretical model, we need to flip the sign of the integral + + # Sort by y first, then z and z (group by y-levels) + sort_idx_top = np.lexsort( + (coords_top[:, 0], -coords_top[:, 1]) + ) # y is sorted, from large to small + sort_idx_bot = np.lexsort( + (coords_bot[:, 0], -coords_bot[:, 1]) + ) # y is sorted, from large to small + + top_coords_sorted = coords_top[ + sort_idx_top + ] # sorted from largest y to smallest y + bot_coords_sorted = coords_bot[ + sort_idx_bot + ] # sorted from largest y to smallest y + + ### integral of top surface: + # Unique x-values (up to numerical tolerance) + ys_top = np.unique(np.round(top_coords_sorted[:, 1], decimals=6)) + + # get f(x) = integral of torque on top surface from left end to x (x is all unique x coordinates on top surface) + indicator_top = 0 + for y_value in ys_top: + # Mask for points at this y-level + y_mask_top = np.isclose(top_coords_sorted[:, 1], y_value, atol=1e-6) + + facet_centers_top = dolfinx.mesh.compute_midpoints( + whole_top_surf_submesh, + 2, + np.array( + np.arange( + whole_top_surf_submesh.topology.index_map(2).size_local + ), + dtype=np.int32, + ), ) - self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid( - domain.domain_markers[f"top_{panel_id:.0f}"]["idx"] + # Identify cells where the midpoint y-coordinate is higher + marked_facet_top = np.where(facet_centers_top[:, 1] >= y_value)[0] + + # Create a submesh from those cells + submesh_top_surface_part, entity_map, vertex_map, geom_map = ( + dolfinx.mesh.create_submesh( + whole_top_surf_submesh, 2, marked_facet_top.astype(np.int32) + ) ) - self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid( - domain.domain_markers[f"right_{panel_id:.0f}"]["idx"] + top_sub_function_space = dolfinx.fem.FunctionSpace( + submesh_top_surface_part, ("CG", 1) ) - self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid( - domain.domain_markers[f"bottom_{panel_id:.0f}"]["idx"] + top_sub_func = dolfinx.fem.Function(top_sub_function_space) + top_sub_func.interpolate(torque_function_on_panels_top) + + dx_top = ufl.Measure("dx", domain=submesh_top_surface_part) + integrated_top = dolfinx.fem.assemble_scalar( + dolfinx.fem.form(top_sub_func * dx_top) ) - self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid( - domain.domain_markers[f"front_{panel_id:.0f}"]["idx"] + top_integrated_func_along_y[sort_idx_top[y_mask_top]] = integrated_top + indicator_top += 1 + ### integral of bot surface: + # Unique y-values (up to numerical tolerance) + ys_bot = np.unique(np.round(bot_coords_sorted[:, 1], decimals=6)) + indicator_bot = 0 + for y_value in ys_bot: + + # Mask for points at this y-level + y_mask_bot = np.isclose(bot_coords_sorted[:, 1], y_value, atol=1e-6) + facet_centers_bot = dolfinx.mesh.compute_midpoints( + whole_bot_surf_submesh, + 2, + np.array( + np.arange( + whole_bot_surf_submesh.topology.index_map(2).size_local + ), + dtype=np.int32, + ), ) - self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid( - domain.domain_markers[f"back_{panel_id:.0f}"]["idx"] + # Identify cells where the midpoint y-coordinate is less than 1 + marked_facet_bot = np.where(facet_centers_bot[:, 1] >= y_value)[0] + # Create a submesh from those cells + submesh_bot_surface_part, entity_map, vertex_map, geom_map = ( + dolfinx.mesh.create_submesh( + whole_bot_surf_submesh, 2, marked_facet_bot.astype(np.int32) + ) + ) + bot_sub_function_space = dolfinx.fem.FunctionSpace( + submesh_bot_surface_part, ("CG", 1) + ) + bot_sub_func = dolfinx.fem.Function(bot_sub_function_space) + bot_sub_func.interpolate(torque_function_on_panels_bot) + + dx_bot = ufl.Measure("dx", domain=submesh_bot_surface_part) + integrated_bot = dolfinx.fem.assemble_scalar( + dolfinx.fem.form(bot_sub_func * dx_bot) + ) + bot_integrated_func_along_y[sort_idx_bot[y_mask_bot]] = integrated_bot + + indicator_bot += 1 + single_integral_function_top = dolfinx.fem.Function( + whole_top_submesh_function_space, name="single_integral_function_top" + ) + single_integral_function_top.x.array[:] = top_integrated_func_along_y + + single_integral_function_bot = dolfinx.fem.Function( + whole_bot_submesh_function_space, name="single_integral_function_bot" + ) + single_integral_function_bot.x.array[:] = bot_integrated_func_along_y + + torque_double_integral = 0.0 + + attr_name = f"double_integral_total_torque_panel_{panel_id:.0f}_0" + # setattr(self, attr_name, torque_double_integral) + + double_integral_total_torque_panel_constant = getattr(self, attr_name) + double_integral_total_torque_panel_constant.value = torque_double_integral + + single_integral_function_top_fluid = dolfinx.fem.Function(self.Q) + single_integral_function_top_fluid.interpolate(single_integral_function_top) + single_integral_function_bot_fluid = dolfinx.fem.Function(self.Q) + single_integral_function_bot_fluid.interpolate(single_integral_function_bot) + + for connector_id in range(params.pv_array.modules_per_span): + torque_double_integral += ( + dolfinx.fem.assemble_scalar( + dolfinx.fem.form( + single_integral_function_top_fluid + * ds_fluid( + domain.domain_markers[ + f"panel_top_{panel_id:.0f}_{params.pv_array.modules_per_span-1-connector_id:.0f}" + ]["idx"] + ) + ) + ) + / params.pv_array.panel_chord + ) + torque_double_integral += ( + dolfinx.fem.assemble_scalar( + dolfinx.fem.form( + single_integral_function_bot_fluid + * ds_fluid( + domain.domain_markers[ + f"panel_bottom_{panel_id:.0f}_{params.pv_array.modules_per_span-1-connector_id:.0f}" + ]["idx"] + ) + ) + ) + / params.pv_array.panel_chord ) + attr_name = f"double_integral_total_torque_panel_{panel_id:.0f}_{connector_id+1:.0f}" - self.integrated_force_z_form[-1] = dolfinx.fem.form( - self.integrated_force_z_form[-1] + double_integral_total_torque_panel_constant = getattr(self, attr_name) + double_integral_total_torque_panel_constant.value = ( + torque_double_integral ) + # double_integral_total_torque_panel_0_0 is the double integral of first panel array at the most back connector (maximum y), + # double_integral_total_torque_panel_0_10 is the double integral of first panel array at the most front connector (smallest y) + def _assemble_system(self, params): """Pre-assemble all LHS matrices and RHS vectors @@ -750,9 +1204,10 @@ def solve(self, domain, params, current_time): self.mesh_vel.vector.array[:] = self.mesh_vel.vector.array / params.solver.dt self.mesh_vel.x.scatter_forward() + # self.compute_double_integral_panel_torques(domain, params) + # Calculate the tentative velocity self._solver_step_1(params) - # Calculate the change in pressure to enforce incompressibility self._solver_step_2(params) @@ -770,6 +1225,10 @@ def solve(self, domain, params, current_time): self.compute_lift_and_drag(params, current_time) + # self.compute_panel_torques(domain, params) + if domain.modeling_torque_tube and params.general.geometry_module == "panels3d": + self.compute_double_integral_panel_torques(domain, params) + # Compute the pressure drop between the inlet and outlet if params.pv_array.stream_rows > 0: self.compute_pressure_drop_between_points(domain, params) diff --git a/pvade/fluid/boundary_conditions.py b/pvade/fluid/boundary_conditions.py index a68b1533..3e6f5dd9 100644 --- a/pvade/fluid/boundary_conditions.py +++ b/pvade/fluid/boundary_conditions.py @@ -585,27 +585,27 @@ def __call__(self, x): or params.general.geometry_module == "heliostats3d" or params.general.geometry_module == "flag2d" ): + for module_id in range(params.pv_array.modules_per_span): + for location in ( + f"panel_bottom_{panel_id}_{module_id}", + f"panel_top_{panel_id}_{module_id}", + f"panel_left_{panel_id}", + f"panel_right_{panel_id}", + # f"front_{panel_id}", # not valid in panels2d? + # f"back_{panel_id}", + ): + T0_pv_panel_scalar = dolfinx.fem.Constant( + domain.fluid.msh, PETSc.ScalarType(params.fluid.T0_panel) + ) - for location in ( - f"bottom_{panel_id}", - f"top_{panel_id}", - f"left_{panel_id}", - f"right_{panel_id}", - # f"front_{panel_id}", # not valid in panels2d? - # f"back_{panel_id}", - ): - T0_pv_panel_scalar = dolfinx.fem.Constant( - domain.fluid.msh, PETSc.ScalarType(params.fluid.T0_panel) - ) - - panel_sfc_dofs = get_facet_dofs_by_gmsh_tag( - domain, functionspace, location - ) - bc = dolfinx.fem.dirichletbc( - T0_pv_panel_scalar, panel_sfc_dofs, functionspace - ) + panel_sfc_dofs = get_facet_dofs_by_gmsh_tag( + domain, functionspace, location + ) + bc = dolfinx.fem.dirichletbc( + T0_pv_panel_scalar, panel_sfc_dofs, functionspace + ) - bcT.append(bc) + bcT.append(bc) # if params.general.debug_flag == True: # print('built temperature boundary conditions') diff --git a/pvade/geometry/MeshManager.py b/pvade/geometry/MeshManager.py index 13f1764d..15408075 100644 --- a/pvade/geometry/MeshManager.py +++ b/pvade/geometry/MeshManager.py @@ -129,6 +129,8 @@ def _get_domain_markers(self, params): def build(self, params): """This function call builds the geometry, marks the boundaries and creates a mesh using Gmsh.""" + self.modeling_torque_tube = False + domain_creation_module = ( f"pvade.geometry.{params.general.geometry_module}.DomainCreation" ) @@ -151,6 +153,8 @@ def build(self, params): and params.general.structural_analysis == True ): self.geometry.build_FSI(params) + self.modeling_torque_tube = self.geometry.modeling_torque_tube + elif ( ( params.general.geometry_module == "panels3d" @@ -160,8 +164,12 @@ def build(self, params): and params.general.structural_analysis == True ): self.geometry.build_structure(params) + self.modeling_torque_tube = self.geometry.modeling_torque_tube else: self.geometry.build_FSI(params) + + self.modeling_torque_tube = False + # Build the domain markers for each surface and cell if hasattr(self.geometry, "domain_markers"): # If the "build" process created domain markers, use those directly... @@ -187,6 +195,8 @@ def build(self, params): self.geometry.ndim ) # gmsh_model.get_dimension() # ?? should this be domain.ndim? + self.modeling_torque_tube = self.comm.bcast(self.modeling_torque_tube, root=0) + # When finished, rank 0 needs to tell other ranks about how the domain_markers dictionary was created # and what values it holds. This is important now since the number of indices "idx" generated in the # geometry module differs from what's prescribed in the init of this class. @@ -214,6 +224,11 @@ def build(self, params): self.cell_tags.name = "cell_tags" self.facet_tags.name = "facet_tags" + with dolfinx.io.XDMFFile(self.comm, f"parent_mesh_only.xdmf", "w") as fp: + fp.write_mesh(self.msh) + fp.write_meshtags(self.cell_tags) + fp.write_meshtags(self.facet_tags) + # if ( # params.general.geometry_module == "panels3d" # or params.general.geometry_module == "heliostats3d" @@ -276,10 +291,30 @@ def _create_submeshes_from_parent(self, params): print(f"Creating {sub_domain_name} submesh") # Get the idx associated with either "fluid" or "structure" - marker_id = self.domain_markers[sub_domain_name]["idx"] - # Find all cells where cell tag = marker_id - submesh_cells = self.cell_tags.find(marker_id) + # if structure includes modules and connectors + if ( + sub_domain_name == "structure" + and "structure" not in self.domain_markers + ): + marker_id = self.domain_markers["modules"]["idx"] + # Find all cells where cell tag = marker_id + if ( + self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + ): + submesh_cells_modules = self.cell_tags.find(marker_id) + marker_id = self.domain_markers["connectors"]["idx"] + submesh_cells = np.hstack( + (self.cell_tags.find(marker_id), submesh_cells_modules) + ) + else: + submesh_cells = self.cell_tags.find(marker_id) + + else: + marker_id = self.domain_markers[sub_domain_name]["idx"] + # Find all cells where cell tag = marker_id + submesh_cells = self.cell_tags.find(marker_id) # Use those found cells to construct a new mesh submesh, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh( @@ -309,7 +344,7 @@ def _transfer_mesh_tags_to_submeshes(self, params): num_facets = f_map.size_local + f_map.num_ghosts all_values = np.zeros(num_facets, dtype=np.int32) - # Assign non-zero facet tags using the facet tag indices + # Assign non-zero facet tags using the facet tag indices, save the facet marker to all_values all_values[self.facet_tags.indices] = self.facet_tags.values cell_to_facet = self.msh.topology.connectivity(self.ndim, facet_dim) @@ -345,14 +380,29 @@ def _transfer_mesh_tags_to_submeshes(self, params): for child, parent in zip(child_facets, parent_facets): sub_values[child] = all_values[parent] + # sub_cell_map = sub_domain.msh.topology.index_map(self.ndim) + f_map_cell = self.msh.topology.index_map(self.ndim) + + # Get the total number of cells in the parent mesh + num_cells = f_map_cell.size_local + f_map_cell.num_ghosts + all_cell_values = np.zeros(num_cells, dtype=np.int32) + all_cell_values[self.cell_tags.indices] = self.cell_tags.values + sub_cell_map = sub_domain.msh.topology.index_map(self.ndim) sub_num_cells = sub_cell_map.size_local + sub_cell_map.num_ghosts + sub_cell_values = np.empty(sub_num_cells, dtype=np.int32) + + for k, entity in enumerate(sub_domain.entity_map): + sub_cell_values[k] = all_cell_values[entity] + + # sub_num_cells = sub_cell_map.size_local + sub_cell_map.num_ghosts + sub_domain.cell_tags = dolfinx.mesh.meshtags( sub_domain.msh, sub_domain.msh.topology.dim, np.arange(sub_num_cells, dtype=np.int32), - np.ones(sub_num_cells, dtype=np.int32), + sub_cell_values, ) sub_domain.cell_tags.name = "cell_tags" @@ -555,6 +605,17 @@ def read_mesh_files(self, read_mesh_dir, params): and use it to solve the CFD/CSD problem """ + if ( + params.pv_array.torque_tube_separation > 0.0 + and params.pv_array.torque_tube_outer_radius > 0.0 + ): + self.modeling_torque_tube = True + else: + self.modeling_torque_tube = False + assert ( + params.pv_array.modules_per_span == 1 + ), "When not modeling torque tube, modules_per_span must be 1." + sub_domain_list = ["fluid", "structure"] for sub_domain_name in sub_domain_list: @@ -713,7 +774,7 @@ def test_mesh_functionspace(self): print(f"Rank {self.rank} owns {num_nodes_owned_by_proc} nodes\n{coords}") - def test_submesh_transfer(self, params): + def test_submesh_transfer(self, params, domain, elasticity): P2 = ufl.VectorElement("Lagrange", self.msh.ufl_cell(), 2) # P2 = ufl.FiniteElement("Lagrange", self.fluid.msh.ufl_cell(), 1) diff --git a/pvade/geometry/flag2d/DomainCreation.py b/pvade/geometry/flag2d/DomainCreation.py index 40eef9ff..352ac0bf 100644 --- a/pvade/geometry/flag2d/DomainCreation.py +++ b/pvade/geometry/flag2d/DomainCreation.py @@ -116,10 +116,10 @@ def build_FSI(self, params): elif np.allclose(com[1], params.domain.y_max): self._add_to_domain_markers("y_max", [surf_id], "facet") - self._add_to_domain_markers("left_0", [5], "facet") - self._add_to_domain_markers("bottom_0", [6], "facet") - self._add_to_domain_markers("right_0", [7], "facet") - self._add_to_domain_markers("top_0", [8], "facet") + self._add_to_domain_markers("panel_left_0", [5], "facet") + self._add_to_domain_markers("panel_bottom_0_0", [6], "facet") + self._add_to_domain_markers("panel_right_0", [7], "facet") + self._add_to_domain_markers("panel_top_0_0", [8], "facet") # Tag objects as either structure or fluid vol_tag_list = self.gmsh_model.occ.getEntities(self.ndim) diff --git a/pvade/geometry/heliostats3d/DomainCreation.py b/pvade/geometry/heliostats3d/DomainCreation.py index c7a804e8..3063d236 100644 --- a/pvade/geometry/heliostats3d/DomainCreation.py +++ b/pvade/geometry/heliostats3d/DomainCreation.py @@ -83,6 +83,8 @@ def Rz(theta): return rot_matrix + self.modeling_torque_tube = False + # Compute and store some useful geometric quantities self.x_span = params.domain.x_max - params.domain.x_min self.y_span = params.domain.y_max - params.domain.y_min @@ -247,22 +249,22 @@ def Rz(theta): # self._add_to_domain_markers(f"z_max_{panel_ct:.0f}", [panel_surfs[-1]], "facet") self._add_to_domain_markers( - f"front_{panel_ct:.0f}", [panel_surfs[0]], "facet" + f"panel_front_{panel_ct:.0f}", [panel_surfs[0]], "facet" ) self._add_to_domain_markers( - f"back_{panel_ct:.0f}", [panel_surfs[1]], "facet" + f"panel_back_{panel_ct:.0f}", [panel_surfs[1]], "facet" ) self._add_to_domain_markers( - f"left_{panel_ct:.0f}", [panel_surfs[2]], "facet" + f"panel_left_{panel_ct:.0f}", [panel_surfs[2]], "facet" ) self._add_to_domain_markers( - f"right_{panel_ct:.0f}", [panel_surfs[3]], "facet" + f"panel_right_{panel_ct:.0f}", [panel_surfs[3]], "facet" ) self._add_to_domain_markers( - f"bottom_{panel_ct:.0f}", panel_surfs[4:-1], "facet" + f"panel_bottom_{panel_ct:.0f}_0", panel_surfs[4:-1], "facet" ) self._add_to_domain_markers( - f"top_{panel_ct:.0f}", [panel_surfs[-1]], "facet" + f"panel_top_{panel_ct:.0f}_0", [panel_surfs[-1]], "facet" ) # self._add_to_domain_markers(f"right_{panel_ct:.0f}", [panel_surfs[1]], "facet")#correct @@ -990,22 +992,22 @@ def set_length_scales(self, params, domain_markers): for panel_id in range(params.pv_array.stream_rows * params.pv_array.span_rows): internal_surface_tags.append( - domain_markers[f"bottom_{panel_id}"]["gmsh_tags"][0] + domain_markers[f"panel_bottom_{panel_id}_0"]["gmsh_tags"][0] ) internal_surface_tags.append( - domain_markers[f"top_{panel_id}"]["gmsh_tags"][0] + domain_markers[f"panel_top_{panel_id}_0"]["gmsh_tags"][0] ) internal_surface_tags.append( - domain_markers[f"left_{panel_id}"]["gmsh_tags"][0] + domain_markers[f"panel_left_{panel_id}"]["gmsh_tags"][0] ) internal_surface_tags.append( - domain_markers[f"right_{panel_id}"]["gmsh_tags"][0] + domain_markers[f"panel_right_{panel_id}"]["gmsh_tags"][0] ) internal_surface_tags.append( - domain_markers[f"front_{panel_id}"]["gmsh_tags"][0] + domain_markers[f"panel_front_{panel_id}"]["gmsh_tags"][0] ) internal_surface_tags.append( - domain_markers[f"back_{panel_id}"]["gmsh_tags"][0] + domain_markers[f"panel_back_{panel_id}"]["gmsh_tags"][0] ) min_dist = [] diff --git a/pvade/geometry/panels2d/DomainCreation.py b/pvade/geometry/panels2d/DomainCreation.py index 0df35acb..9c4a30b9 100644 --- a/pvade/geometry/panels2d/DomainCreation.py +++ b/pvade/geometry/panels2d/DomainCreation.py @@ -184,22 +184,22 @@ def build_FSI(self, params): if np.isclose(com[0], x_min_panel): self._add_to_domain_markers( - f"left_{panel_id:.0f}", [surf_id], "facet" + f"panel_left_{panel_id:.0f}", [surf_id], "facet" ) elif np.isclose(com[0], x_max_panel): self._add_to_domain_markers( - f"right_{panel_id:.0f}", [surf_id], "facet" + f"panel_right_{panel_id:.0f}", [surf_id], "facet" ) elif np.isclose(com[1], y_min_panel): self._add_to_domain_markers( - f"bottom_{panel_id:.0f}", [100], "facet" + f"panel_bottom_{panel_id:.0f}_0", [surf_id], "facet" ) elif np.isclose(com[1], y_max_panel): self._add_to_domain_markers( - f"top_{panel_id:.0f}", [surf_id], "facet" + f"panel_top_{panel_id:.0f}_0", [surf_id], "facet" ) # Rotate the panel currently centered at diff --git a/pvade/geometry/panels3d/DomainCreation.py b/pvade/geometry/panels3d/DomainCreation.py index 5ead1e4b..540870b3 100644 --- a/pvade/geometry/panels3d/DomainCreation.py +++ b/pvade/geometry/panels3d/DomainCreation.py @@ -112,6 +112,17 @@ def Rz(theta): array_rotation = (params.fluid.wind_direction + 90.0) % 360.0 array_rotation_rad = np.radians(array_rotation) + if ( + params.pv_array.torque_tube_separation > 0.0 + and params.pv_array.torque_tube_outer_radius > 0.0 + ): + self.modeling_torque_tube = True + else: + self.modeling_torque_tube = False + assert ( + params.pv_array.modules_per_span == 1 + ), "When not modeling torque tube, modules_per_span must be 1." + # The centroid of each panel in the x-direction (these should start at x=0) x_centers = np.linspace( 0.0, @@ -149,11 +160,25 @@ def Rz(theta): domain_tag_list = [] domain_tag_list.append(domain_tag) + # panel_tag_list includes all structure, panel+connector panel_tag_list = [] panel_ct = 0 - prev_surf_tag = [] + # only include panels + structure_panel_only_list = [] + # only include connectors + structure_connector_only_list = [] + + # if params.pv_array.modules_per_span > 1: + module_distances = np.linspace( + -half_span, half_span, params.pv_array.modules_per_span + 1 + ) + + transformed_com = {} # all panels + vol_tags_modules = [] + vol_tags_connectors = [] + # start to add panel for panel_id_y, yy in enumerate(y_centers): for panel_id_x, xx in enumerate(x_centers): if isinstance(params.pv_array.tracker_angle, list): @@ -169,205 +194,797 @@ def Rz(theta): else: tracker_angle_rad = np.radians(params.pv_array.tracker_angle) - # Create an 0-tracking-degree panel centered at (x, y, z) = (0, 0, 0) - panel_id = self.gmsh_model.occ.addBox( - -half_chord, - -half_span, - -half_thickness, - params.pv_array.panel_chord, - params.pv_array.panel_span, - params.pv_array.panel_thickness, - ) + this_panel_tag_list = [] # for each row + this_panel_transformed_com = {} # for each row - panel_tag = (self.ndim, panel_id) - panel_tag_list.append(panel_tag) + numpy_pt_list = [] # for each row + embedded_lines_tag_list = [] # for each row - numpy_pt_list = [] - embedded_lines_tag_list = [] + for module_id in range(params.pv_array.modules_per_span): - # Add a bisecting line to the bottom of the panel in the spanwise direction - pt_1 = self.gmsh_model.occ.addPoint(0, -half_span, -half_thickness) - pt_2 = self.gmsh_model.occ.addPoint(0, half_span, -half_thickness) - - numpy_pt_list.append( - [0, -half_span, -half_thickness, 0, half_span, -half_thickness] - ) + module_span = ( + module_distances[module_id + 1] - module_distances[module_id] + ) - torque_tube_id = self.gmsh_model.occ.addLine(pt_1, pt_2) - torque_tube_tag = (1, torque_tube_id) - embedded_lines_tag_list.append(torque_tube_tag) - - # Add lines in the streamwise direction to mimic sections of panel held rigid by motor - if params.pv_array.span_fixation_pts is not None: - if not isinstance(params.pv_array.span_fixation_pts, list): - num_fixation_pts = int( - np.floor( - params.pv_array.panel_span - / params.pv_array.span_fixation_pts - ) + if ( + self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + ): + # Create an 0-tracking-degree panel centered at (x, y, z) = (0, 0, 0) + this_module = self.gmsh_model.occ.addBox( + -half_chord, + module_distances[module_id], + params.pv_array.torque_tube_separation - half_thickness, + params.pv_array.panel_chord, + module_span, + params.pv_array.panel_thickness, ) - fixation_pts_list = [] - - for k in range(1, num_fixation_pts + 1): - next_pt = k * params.pv_array.span_fixation_pts - - eps = 1e-9 - - if ( - next_pt > eps - and next_pt < params.pv_array.panel_span - eps - ): - fixation_pts_list.append(next_pt) + panel_tag_list.append((self.ndim, this_module)) # all panels + this_panel_tag_list.append( + (self.ndim, this_module) + ) # this panel array/row + + structure_panel_only_list.append((self.ndim, this_module)) + + this_standoff = self.gmsh_model.occ.addBox( + -params.pv_array.block_chord_div_by_panel_chord + * half_chord, + module_distances[module_id], + -half_thickness, + 2.0 + * params.pv_array.block_chord_div_by_panel_chord + * half_chord, + params.pv_array.block_chord_div_by_panel_chord * half_chord, + params.pv_array.torque_tube_separation, + ) + panel_tag_list.append((self.ndim, this_standoff)) + this_panel_tag_list.append((self.ndim, this_standoff)) - else: - fixation_pts_list = params.pv_array.span_fixation_pts + structure_connector_only_list.append((self.ndim, this_standoff)) - for fp in fixation_pts_list: + # Add a bisecting line to the bottom of the connector in the spanwise direction pt_1 = self.gmsh_model.occ.addPoint( - -half_chord, -half_span + fp, -half_thickness + 0, module_distances[module_id], -half_thickness ) pt_2 = self.gmsh_model.occ.addPoint( - half_chord, -half_span + fp, -half_thickness + 0, + module_distances[module_id] + + params.pv_array.block_chord_div_by_panel_chord + * half_chord, + -half_thickness, ) - - # FIXME: don't add the fixation points into the numpy tagging for now numpy_pt_list.append( [ - -half_chord, - -half_span + fp, + 0, + module_distances[module_id], -half_thickness, - half_chord, - -half_span + fp, + 0, + module_distances[module_id] + + params.pv_array.block_chord_div_by_panel_chord + * half_chord, -half_thickness, ] + ) # for this row + torque_tube_id = self.gmsh_model.occ.addLine(pt_1, pt_2) + torque_tube_tag = (1, torque_tube_id) + embedded_lines_tag_list.append( + torque_tube_tag + ) # for this panel row + + else: # this loop is to add panel only, no torque tube modeling + this_module = self.gmsh_model.occ.addBox( + -half_chord, + module_distances[module_id], + -half_thickness, + params.pv_array.panel_chord, + module_span, + params.pv_array.panel_thickness, ) + panel_tag_list.append((self.ndim, this_module)) + this_panel_tag_list.append((self.ndim, this_module)) - fixed_pt_id = self.gmsh_model.occ.addLine(pt_1, pt_2) - fixed_pt_tag = (1, fixed_pt_id) - - embedded_lines_tag_list.append(fixed_pt_tag) + structure_panel_only_list.append((self.ndim, this_module)) - # Store the result of fragmentation, it holds all the small surfaces we need to tag - panel_frags = self.gmsh_model.occ.fragment( - [panel_tag], embedded_lines_tag_list - ) - - # extract just the first entry, and remove the 3d entry in position 0 - panel_surfs = panel_frags[0] - panel_surfs.pop(0) - panel_surfs = [k[1] for k in panel_surfs] - - # TODO: USE THESE UNAMBIGUOUS NAMES IN A FUTURE REFACTOR - # self._add_to_domain_markers(f"x_min_{panel_ct:.0f}", [panel_surfs[0]], "facet") - # self._add_to_domain_markers(f"x_max_{panel_ct:.0f}", [panel_surfs[1]], "facet") - # self._add_to_domain_markers(f"y_min_{panel_ct:.0f}", [panel_surfs[2]], "facet") - # self._add_to_domain_markers(f"y_max_{panel_ct:.0f}", [panel_surfs[3]], "facet") - # self._add_to_domain_markers(f"z_min_{panel_ct:.0f}", panel_surfs[4:-1], "facet") - # self._add_to_domain_markers(f"z_max_{panel_ct:.0f}", [panel_surfs[-1]], "facet") + # Add a bisecting line to the bottom of the panel in the spanwise direction + pt_1 = self.gmsh_model.occ.addPoint( + 0, module_distances[module_id], -half_thickness + ) + pt_2 = self.gmsh_model.occ.addPoint( + 0, module_distances[module_id + 1], -half_thickness + ) + # numpy_pt_list for this row + numpy_pt_list.append( + [ + 0, + module_distances[module_id], + -half_thickness, + 0, + module_distances[module_id + 1], + -half_thickness, + ] + ) + torque_tube_id = self.gmsh_model.occ.addLine( + pt_1, pt_2 + ) # line simulating torque tube attached horizontally at the bottom of the module + torque_tube_tag = (1, torque_tube_id) + # Store all the embedded lines that we need in the mesh + embedded_lines_tag_list.append( + torque_tube_tag + ) # for this panel row + + # Add lines in the streamwise direction to mimic sections of panel held rigid by motor + if params.pv_array.span_fixation_pts is not None: + if not isinstance(params.pv_array.span_fixation_pts, list): + num_fixation_pts = int( + np.floor( + params.pv_array.panel_span + / params.pv_array.span_fixation_pts + ) + ) + + fixation_pts_list = [] + + for k in range(1, num_fixation_pts + 1): + next_pt = k * params.pv_array.span_fixation_pts + + eps = 1e-9 + + if ( + next_pt > eps + and next_pt < params.pv_array.panel_span - eps + ): + + fixation_pts_list.append(next_pt) + + else: + fixation_pts_list = params.pv_array.span_fixation_pts + + for fp in fixation_pts_list: + pt_1 = self.gmsh_model.occ.addPoint( + -half_chord, -half_span + fp, -half_thickness + ) + pt_2 = self.gmsh_model.occ.addPoint( + half_chord, -half_span + fp, -half_thickness + ) + + # FIXME: don't add the fixation points into the numpy tagging for now + numpy_pt_list.append( + [ + -half_chord, + -half_span + fp, + -half_thickness, + half_chord, + -half_span + fp, + -half_thickness, + ] + ) + + fixed_pt_id = self.gmsh_model.occ.addLine(pt_1, pt_2) + fixed_pt_tag = (1, fixed_pt_id) + + embedded_lines_tag_list.append(fixed_pt_tag) + + if ( + self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + ): # add the last connector + last_standoff = self.gmsh_model.occ.addBox( + -params.pv_array.block_chord_div_by_panel_chord * half_chord, + module_distances[module_id + 1] + - params.pv_array.block_chord_div_by_panel_chord * half_chord, + -half_thickness, + 2.0 + * params.pv_array.block_chord_div_by_panel_chord + * half_chord, + params.pv_array.block_chord_div_by_panel_chord * half_chord, + params.pv_array.torque_tube_separation, + ) + this_panel_tag_list.append((self.ndim, last_standoff)) + panel_tag_list.append((self.ndim, last_standoff)) - # Translate the panel by (x_center, y_center, elev) - self.gmsh_model.occ.translate( - [panel_tag], - xx, - yy, - params.pv_array.elevation, - ) + structure_connector_only_list.append((self.ndim, last_standoff)) - # Get the bounding box for this panel - panel_bounding_box = self.gmsh_model.occ.getBoundingBox( - panel_tag[0], panel_tag[1] + pt_1 = self.gmsh_model.occ.addPoint( + 0, + module_distances[module_id + 1] + - params.pv_array.block_chord_div_by_panel_chord * half_chord, + -half_thickness, + ) + pt_2 = self.gmsh_model.occ.addPoint( + 0, module_distances[module_id + 1], -half_thickness + ) + numpy_pt_list.append( + [ + 0, + module_distances[module_id + 1] + - params.pv_array.block_chord_div_by_panel_chord + * half_chord, + -half_thickness, + 0, + module_distances[module_id + 1], + -half_thickness, + ] + ) + torque_tube_id = self.gmsh_model.occ.addLine(pt_1, pt_2) + torque_tube_tag = (1, torque_tube_id) + embedded_lines_tag_list.append(torque_tube_tag) + + """ this is not used as we are using block surface to simulate the fixed bc applied by motor """ + # # Add lines in the streamwise direction to mimic sections of panel held rigid by motor + # if params.pv_array.span_fixation_pts is not None: + # if not isinstance(params.pv_array.span_fixation_pts, list): + # num_fixation_pts = int( + # np.floor( + # params.pv_array.panel_span + # / params.pv_array.span_fixation_pts + # ) + # ) + + # fixation_pts_list = [] + + # for k in range(1, num_fixation_pts + 1): + # next_pt = k * params.pv_array.span_fixation_pts + + # eps = 1e-9 + + # if ( + # next_pt > eps + # and next_pt < params.pv_array.panel_span - eps + # ): + # fixation_pts_list.append(next_pt) + + # else: + # fixation_pts_list = params.pv_array.span_fixation_pts + + # for fp in fixation_pts_list: + # # FIXME: don't add the fixation points into the numpy tagging for now + # if modeling_torque_tube: + # numpy_pt_list.append( + # [ + # -params.pv_array.torque_tube_radius, + # -half_span + fp, + # 0.0, + # params.pv_array.torque_tube_radius, + # -half_span + fp, + # 0.0, + # ] + # ) + + # else: + # numpy_pt_list.append( + # [ + # -half_chord, + # -half_span + fp, + # 0.0, + # half_chord, + # -half_span + fp, + # 0.0, + # ] + # ) + + # # If the separation line already exists at a module division, + # # there's no need to redraw it, but otherwise, add the gmsh points + # # to the model for embedding + # if not np.any(np.isclose(-half_span + fp, module_distances)): + # print( + # f"Embedding a line for a no-deformation boundary condition." + # ) + + # if modeling_torque_tube: + # pt_1 = self.gmsh_model.occ.addPoint( + # -params.pv_array.torque_tube_radius, + # -half_span + fp, + # 0.0, + # ) + # pt_2 = self.gmsh_model.occ.addPoint( + # params.pv_array.torque_tube_radius, + # -half_span + fp, + # 0.0, + # ) + + # else: + # pt_1 = self.gmsh_model.occ.addPoint( + # -half_chord, -half_span + fp, 0.0 + # ) + # pt_2 = self.gmsh_model.occ.addPoint( + # half_chord, -half_span + fp, 0.0 + # ) + + # fixation_line_id = self.gmsh_model.occ.addLine(pt_1, pt_2) + # fixation_line_tag = (1, fixation_line_id) + # embedded_lines_tag_list.append(fixation_line_tag) + # else: + # print( + # f"Applying no-deformation boundary condition at {fp}." + # ) + + # Fragment the lines into the surfaces (equivalent to embedding these lines) + self.gmsh_model.occ.fragment( + this_panel_tag_list, embedded_lines_tag_list ) - # Store the x_min, y_min, x_max, y_max points - # corresponding to the un-rotated (tracking_angle = 0) configuration - x_min_panel = panel_bounding_box[0] - y_min_panel = panel_bounding_box[1] - z_min_panel = panel_bounding_box[2] - x_max_panel = panel_bounding_box[3] - y_max_panel = panel_bounding_box[4] - z_max_panel = panel_bounding_box[5] - - self.gmsh_model.occ.synchronize() + for panel_tag in this_panel_tag_list: # 3d cell domain + self.gmsh_model.occ.synchronize() - # Get the list of 1D surfaces (edges) that make up this panel - surf_tags_for_this_panel = self.gmsh_model.getBoundary( - [panel_tag], oriented=False - ) + # Get the list of 2D surfaces (surfaces) that make up this panel + surf_tags_for_this_panel = self.gmsh_model.getBoundary( + [panel_tag], oriented=False + ) - for surf_tag in surf_tags_for_this_panel: - surf_dim = surf_tag[0] - surf_id = surf_tag[1] - com = self.gmsh_model.occ.getCenterOfMass(surf_dim, surf_id) - # sturctures tagging - if np.isclose(com[0], x_min_panel): - self._add_to_domain_markers( - f"left_{panel_ct:.0f}", [surf_id], "facet" - ) + vol_com = self.gmsh_model.occ.getCenterOfMass( + self.ndim, panel_tag[1] + ) + if np.isclose( + vol_com[2], + params.pv_array.torque_tube_separation, + ): + target_key = f"modules" + elif np.isclose( + vol_com[2], + params.pv_array.torque_tube_separation / 2.0 - half_thickness, + ): + target_key = f"connectors" + + if target_key is not None: + if target_key in this_panel_transformed_com: + this_panel_transformed_com[target_key].append(vol_com) + else: + this_panel_transformed_com[target_key] = [vol_com] + + for surf_tag in surf_tags_for_this_panel: + surf_dim = surf_tag[0] + surf_id = surf_tag[1] + com = self.gmsh_model.occ.getCenterOfMass(surf_dim, surf_id) + + target_key = None + + surface_located_or_not = False + + # sturctures tagging + if np.isclose(com[0], -half_chord): + target_key = f"panel_left_{panel_ct:.0f}" + surface_located_or_not = True + + if np.isclose(com[0], half_chord): + target_key = f"panel_right_{panel_ct:.0f}" + surface_located_or_not = True + + if ( + np.isclose(com[1], module_distances[0]) + and np.isclose(com[0], 0) + and np.isclose( + com[2], + params.pv_array.torque_tube_separation, + ) + ): + target_key = f"panel_front_{panel_ct:.0f}" + surface_located_or_not = True + + if ( + np.isclose( + com[1], + module_distances[params.pv_array.modules_per_span], + ) + and np.isclose(com[0], 0) + and np.isclose( + com[2], + params.pv_array.torque_tube_separation, + ) + ): + target_key = f"panel_back_{panel_ct:.0f}" + surface_located_or_not = True + + if (not self.modeling_torque_tube) or ( + params.general.geometry_module != "panels3d" + ): + if np.isclose( + com[2], + params.pv_array.torque_tube_separation - half_thickness, + ): + target_key = f"panel_bottom_{panel_ct:.0f}_0" + surface_located_or_not = True - elif np.isclose(com[0], x_max_panel): - self._add_to_domain_markers( - f"right_{panel_ct:.0f}", [surf_id], "facet" - ) + if np.isclose( + com[2], + half_thickness + params.pv_array.torque_tube_separation, + ): + target_key = f"panel_top_{panel_ct:.0f}_0" + surface_located_or_not = True + + else: + for module_id in range(params.pv_array.modules_per_span): + if ( + com[1] + >= module_distances[module_id] + + module_span / 2.0 + - params.pv_array.block_chord_div_by_panel_chord + * half_chord + and com[1] + <= module_distances[module_id] + + module_span / 2.0 + + params.pv_array.block_chord_div_by_panel_chord + * half_chord + and np.isclose(com[0], 0) + and np.isclose( + com[2], + params.pv_array.torque_tube_separation + - half_thickness, + ) + ): + target_key = ( + f"panel_bottom_{panel_ct:.0f}_{module_id:.0f}" + ) + surface_located_or_not = True + break + + if ( + com[1] + >= module_distances[module_id] + + module_span / 2.0 + - params.pv_array.block_chord_div_by_panel_chord + * half_chord + and com[1] + <= module_distances[module_id] + + module_span / 2.0 + + params.pv_array.block_chord_div_by_panel_chord + * half_chord + and np.isclose(com[0], 0) + and np.isclose( + com[2], + params.pv_array.torque_tube_separation + + half_thickness, + ) + ): + target_key = ( + f"panel_top_{panel_ct:.0f}_{module_id:.0f}" + ) + surface_located_or_not = True + break + # if ( + # np.isclose( + # com[2], + # params.pv_array.torque_tube_separation + # + params.pv_array.panel_thickness + # ) + + # ): + # target_key = f"panel_top_{panel_ct:.0f}" + # surface_located_or_not = True + + # mark the last block in each row + if ( + np.isclose( + com[0], + params.pv_array.block_chord_div_by_panel_chord + * half_chord, + ) + and self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + and np.isclose( + com[1], + module_distances[params.pv_array.modules_per_span] + - params.pv_array.block_chord_div_by_panel_chord + * half_chord + / 2.0, + ) + ): + target_key = f"block_right_{panel_ct:.0f}_{params.pv_array.modules_per_span:.0f}" + surface_located_or_not = True - elif np.isclose(com[1], y_min_panel): - self._add_to_domain_markers( - f"front_{panel_ct:.0f}", [surf_id], "facet" - ) + if ( + np.isclose( + com[0], + -params.pv_array.block_chord_div_by_panel_chord + * half_chord, + ) + and self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + and np.isclose( + com[1], + module_distances[params.pv_array.modules_per_span] + - params.pv_array.block_chord_div_by_panel_chord + * half_chord + / 2.0, + ) + ): + target_key = f"block_left_{panel_ct:.0f}_{params.pv_array.modules_per_span:.0f}" + surface_located_or_not = True - elif np.isclose(com[1], y_max_panel): - self._add_to_domain_markers( - f"back_{panel_ct:.0f}", [surf_id], "facet" - ) + if ( + np.isclose( + com[2], + params.pv_array.torque_tube_separation / 2.0 + - half_thickness, + ) + and self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + and np.isclose( + com[1], + module_distances[params.pv_array.modules_per_span] + - params.pv_array.block_chord_div_by_panel_chord + * half_chord, + ) + ): + target_key = f"block_front_{panel_ct:.0f}_{params.pv_array.modules_per_span:.0f}" + surface_located_or_not = True - elif np.isclose(com[2], z_min_panel): - self._add_to_domain_markers( - f"bottom_{panel_ct:.0f}", [surf_id], "facet" - ) + if ( + np.isclose( + com[2], + params.pv_array.torque_tube_separation / 2.0 + - half_thickness, + ) + and self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + and np.isclose( + com[1], + module_distances[params.pv_array.modules_per_span], + ) + ): + target_key = f"block_back_{panel_ct:.0f}_{params.pv_array.modules_per_span:.0f}" + surface_located_or_not = True - elif np.isclose(com[2], z_max_panel): - self._add_to_domain_markers( - f"top_{panel_ct:.0f}", [surf_id], "facet" - ) + if ( + np.isclose(com[2], -half_thickness) + and self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + and np.isclose( + com[1], + module_distances[params.pv_array.modules_per_span] + - params.pv_array.block_chord_div_by_panel_chord + * half_chord + / 2.0, + ) + ): + target_key = f"block_bottom_{panel_ct:.0f}_{params.pv_array.modules_per_span:.0f}" + surface_located_or_not = True + if ( + np.isclose( + com[2], + params.pv_array.torque_tube_separation + - half_thickness, + ) + and self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + and np.isclose( + com[1], + module_distances[params.pv_array.modules_per_span] + - params.pv_array.block_chord_div_by_panel_chord + * half_chord + / 2.0, + ) + ): + target_key = f"interior_surface_{panel_ct:.0f}" # block/panel interface + surface_located_or_not = True + + # if not the most left/right panel boundary, it is panel/panel interface + if not surface_located_or_not: + for module_id in range( + params.pv_array.modules_per_span + ): + + # sturctures tagging + + if ( + np.isclose(com[1], module_distances[module_id]) + and np.isclose(com[0], 0) + and np.isclose( + com[2], + params.pv_array.torque_tube_separation, + ) + ): + target_key = f"interior_surface_{panel_ct:.0f}" + surface_located_or_not = True + break + + if ( + np.isclose( + com[1], module_distances[module_id + 1] + ) + and np.isclose(com[0], 0) + and np.isclose( + com[2], + params.pv_array.torque_tube_separation, + ) + ): + target_key = f"interior_surface_{panel_ct:.0f}" + surface_located_or_not = True + break + + # if ( + # np.isclose(com[2], params.pv_array.torque_tube_separation) + # and np.isclose(com[0], 0.0) and com[1] >= module_distances[module_id]+module_span/2.0-params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0 + # and com[1] <= module_distances[module_id]+module_span/2.0+params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0 + # ): + # target_key = f"panel_bottom_{panel_ct:.0f}" + # surface_located_or_not = True + # break + + if ( + np.isclose( + com[2], + params.pv_array.torque_tube_separation + - half_thickness, + ) + and np.isclose(com[0], 0.0) + and np.isclose( + com[1], + module_distances[module_id] + + params.pv_array.block_chord_div_by_panel_chord + * half_chord + / 2.0, + ) + and self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + ): + target_key = f"interior_surface_{panel_ct:.0f}" + surface_located_or_not = True + break + + if ( + np.isclose( + com[0], + params.pv_array.block_chord_div_by_panel_chord + * half_chord, + ) + and self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + and np.isclose( + com[1], + module_distances[module_id] + + params.pv_array.block_chord_div_by_panel_chord + * half_chord + / 2.0, + ) + ): + target_key = f"block_right_{panel_ct:.0f}_{module_id:.0f}" + surface_located_or_not = True + break + + if ( + np.isclose( + com[0], + -params.pv_array.block_chord_div_by_panel_chord + * half_chord, + ) + and self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + and np.isclose( + com[1], + module_distances[module_id] + + params.pv_array.block_chord_div_by_panel_chord + * half_chord + / 2.0, + ) + ): + target_key = ( + f"block_left_{panel_ct:.0f}_{module_id:.0f}" + ) + surface_located_or_not = True + break + + if ( + np.isclose( + com[2], + params.pv_array.torque_tube_separation / 2.0 + - half_thickness, + ) + and self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + and np.isclose( + com[1], module_distances[module_id] + ) + ): + target_key = f"block_front_{panel_ct:.0f}_{module_id:.0f}" + surface_located_or_not = True + break + + if ( + np.isclose( + com[2], + params.pv_array.torque_tube_separation / 2.0 + - half_thickness, + ) + and self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + and np.isclose( + com[1], + module_distances[module_id] + + params.pv_array.block_chord_div_by_panel_chord + * half_chord, + ) + ): + target_key = ( + f"block_back_{panel_ct:.0f}_{module_id:.0f}" + ) + surface_located_or_not = True + break + + if ( + np.isclose(com[2], -half_thickness) + and self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + and np.isclose( + com[1], + module_distances[module_id] + + params.pv_array.block_chord_div_by_panel_chord + * half_chord + / 2.0, + ) + ): + target_key = f"block_bottom_{panel_ct:.0f}_{module_id:.0f}" + surface_located_or_not = True + break + + if not surface_located_or_not: + target_key = f"trash_{panel_ct:.0f}" + print("facet in trash") + + if target_key is not None: + if target_key in this_panel_transformed_com: + this_panel_transformed_com[target_key].append(com) + else: + this_panel_transformed_com[target_key] = [com] + + # print(this_panel_transformed_com[f"trash_{panel_ct:.0f}"]) + + for ( + key, + val, + ) in ( + this_panel_transformed_com.items() + ): # all facets for each panel row. + for row_num, com in enumerate(val): + com_array = np.array(com) + + com_array = np.dot(com_array, Ry(tracker_angle_rad).T) + + com_array[0] += xx + com_array[1] += yy + com_array[2] += params.pv_array.elevation + + com_array[0] -= x_center_of_mass + com_array[1] -= y_center_of_mass + + com_array = np.dot(com_array, Rz(array_rotation_rad).T) + + com_array[0] += x_center_of_mass + com_array[1] += y_center_of_mass + + if key in transformed_com: + transformed_com[key].append(com_array) + else: + transformed_com[key] = [com_array] # including all rows + + # actually this is panel array count panel_ct += 1 # Rotate the panel by its tracking angle along the y-axis - # (currently centered at (xx, yy, params.pv_array.elevation)) + # (currently centered at (0.0, 0.0, 0.0)) self.gmsh_model.occ.rotate( - [panel_tag], - xx, - yy, - params.pv_array.elevation, + this_panel_tag_list, + 0.0, + 0.0, + 0.0, 0, 1, 0, tracker_angle_rad, ) - # Note that the numpy pt panel array has NOT been shifted, i.e., - # it still represents a panel centered at (0, 0, 0), so the rotation - # can be applied directly about the point (0, 0, 0) without any - # shifting like that which needs to be done above - numpy_pt_panel_array = np.array(numpy_pt_list) - numpy_pt_panel_array = np.reshape(numpy_pt_panel_array, (-1, self.ndim)) - - numpy_pt_panel_array = np.dot( - numpy_pt_panel_array, Ry(tracker_angle_rad).T + # Translate the panel by (x_center, y_center, elev) + self.gmsh_model.occ.translate( + this_panel_tag_list, + xx, + yy, + params.pv_array.elevation, ) - # if not hasattr(self, "numpy_pt_array"): - # numpy_pt_array = np.array(numpy_pt_list) - # else: - # numpy_pt_array = np.vcat(numpy_pt_array, np.array(numpy_pt_list)) - - numpy_pt_panel_array[:, 0] += xx - numpy_pt_panel_array[:, 1] += yy - numpy_pt_panel_array[:, 2] += params.pv_array.elevation - # Rotate the panel about the center of the full array as a proxy for changing wind direction (x_center, y_center, 0) self.gmsh_model.occ.rotate( - [panel_tag], + this_panel_tag_list, x_center_of_mass, y_center_of_mass, 0, @@ -377,6 +994,22 @@ def Rz(theta): array_rotation_rad, ) + # Now, apply the same transformations to the numpy representation + # Rotate the panel by its tracking angle along the y-axis + # (currently centered at (0.0, 0.0, 0.0)) + numpy_pt_panel_array = np.array(numpy_pt_list) + numpy_pt_panel_array = np.reshape(numpy_pt_panel_array, (-1, self.ndim)) + + numpy_pt_panel_array = np.dot( + numpy_pt_panel_array, Ry(tracker_angle_rad).T + ) + + # Translate the panel by (x_center, y_center, elev) + numpy_pt_panel_array[:, 0] += xx + numpy_pt_panel_array[:, 1] += yy + numpy_pt_panel_array[:, 2] += params.pv_array.elevation + + # Rotate the panel about the center of the full array as a proxy for changing wind direction (x_center, y_center, 0) numpy_pt_panel_array[:, 0] -= x_center_of_mass numpy_pt_panel_array[:, 1] -= y_center_of_mass @@ -394,49 +1027,30 @@ def Rz(theta): else: self.numpy_pt_total_array = np.copy(numpy_pt_panel_array) - # Check that this panel still exists in the confines of the domain - bbox = self.gmsh_model.occ.get_bounding_box(panel_tag[0], panel_tag[1]) - - if bbox[0] < params.domain.x_min: - raise ValueError( - f"Panel with location (x, y) = ({xx}, {yy}) extends past x_min wall." - ) - if bbox[1] < params.domain.y_min: - raise ValueError( - f"Panel with location (x, y) = ({xx}, {yy}) extends past y_min wall." - ) - if bbox[3] > params.domain.x_max: - raise ValueError( - f"Panel with location (x, y) = ({xx}, {yy}) extends past x_max wall." - ) - if bbox[4] > params.domain.y_max: - raise ValueError( - f"Panel with location (x, y) = ({xx}, {yy}) extends past y_max wall." - ) - # Fragment all panels from the overall domain self.gmsh_model.occ.fragment(domain_tag_list, panel_tag_list) self.gmsh_model.occ.synchronize() + gmsh.write("pnales.brep") + # exit() self.numpy_pt_total_array = np.reshape( self.numpy_pt_total_array, (-1, int(2 * self.ndim)) ) - # import matplotlib.pyplot as plt - # for k in self.numpy_pt_total_array: - # plt.plot([k[0], k[3]], [k[1], k[4]]) - # plt.show() + # print(transformed_com) + # exit() - # it is not necessary to loop over all surfaces, since the panel - # surfaces have been tagged already, but this ensures any ordering change - # doesn't cause problems in the future + # for all panels + # Loop over all the finalized surfaces after fragmentation and tag everything all_surf_tag_list = self.gmsh_model.occ.getEntities(self.ndim - 1) for surf_tag in all_surf_tag_list: surf_id = surf_tag[1] com = self.gmsh_model.occ.getCenterOfMass(self.ndim - 1, surf_id) + located_this_surface = False + # sturctures tagging if np.isclose(com[0], params.domain.x_min): self._add_to_domain_markers("x_min", [surf_id], "facet") @@ -456,22 +1070,91 @@ def Rz(theta): elif np.allclose(com[2], params.domain.z_max): self._add_to_domain_markers("z_max", [surf_id], "facet") + else: + for key, val in transformed_com.items(): + for target_com in val: + # print(target_com) + if np.allclose(np.array(com), target_com): + located_this_surface = True + if "trash" not in key: + # print(key) + self._add_to_domain_markers(key, [surf_id], "facet") + # if "interior_surface" not in key: + # self._add_to_domain_markers("structure_fluid_interface", [surf_id], "facet") + + if not located_this_surface: + print( + f"Warning: Surface {surf_tag} has not been added to domain markers" + ) + + # Since this is not one of the exterior walls, we should check if it extends + # past the boundaries x_min, x_max, ... + this_surf_bbox = self.gmsh_model.occ.get_bounding_box( + self.ndim - 1, surf_id + ) + + # Test that the rotated point still exists in the box domain + if this_surf_bbox[0] < params.domain.x_min: + raise ValueError(f"A panel extends past the x_min wall.") + if this_surf_bbox[0] > params.domain.x_max: + raise ValueError(f"A panel extends past the x_max wall.") + if this_surf_bbox[1] < params.domain.y_min: + raise ValueError(f"A panel extends past the y_min wall.") + if this_surf_bbox[1] > params.domain.y_max: + raise ValueError(f"A panel extends past the y_max wall.") + if this_surf_bbox[2] < 0.0: # params.domain.z_min: + raise ValueError( + f"A panel extends past the z_min wall (ground level = 0.0)." + ) + if this_surf_bbox[2] > params.domain.z_max: + raise ValueError(f"A panel extends past the z_max wall.") + + # mark the panel and connector volumes + all_vol_tag_list = self.gmsh_model.occ.getEntities(self.ndim) + + for vol_tag in all_vol_tag_list: + vol_id = vol_tag[1] + if vol_id != all_vol_tag_list[-1][-1]: + com = self.gmsh_model.occ.getCenterOfMass(self.ndim, vol_id) + + located_this_volume = False + + for key, val in transformed_com.items(): + for target_com in val: + # print(target_com) + if np.allclose(np.array(com), target_com): + located_this_volume = True + if "trash" not in key: + self._add_to_domain_markers(key, [vol_id], "cell") + + # Mark the fluid volume # Volumes are the entities with dimension equal to the mesh dimension vol_tag_list = self.gmsh_model.occ.getEntities(self.ndim) structure_vol_list = [] fluid_vol_list = [] + num_solids = len(vol_tag_list) + for vol_tag in vol_tag_list: vol_id = vol_tag[1] - if vol_id <= params.pv_array.stream_rows * params.pv_array.span_rows: + if vol_id <= num_solids - 1: # Solid Cell + # Since all panel/table structures were created after the box domain + # the final fragmentation removes the original vol_tag of the box and appends + # it to the end, thus, everything with id < num_solids is structure, and + # vol_tag = num_solids (the last-added, fragmented box) is fluid structure_vol_list.append(vol_id) else: # Fluid Cell fluid_vol_list.append(vol_id) - self._add_to_domain_markers("structure", structure_vol_list, "cell") + structure_panel_only_list = [] + structure_connector_only_list = [] + # + # # self._add_to_domain_markers("structure", structure_vol_list, "cell") + # for individual_vol_id in structure_vol_list: + # self._add_to_domain_markers(f"structure_{individual_vol_id:03.0f}", [individual_vol_id], "cell") self._add_to_domain_markers("fluid", fluid_vol_list, "cell") # Record all the data collected to domain_markers as physics groups with physical names @@ -608,6 +1291,17 @@ def Rz(theta): array_rotation = (params.fluid.wind_direction + 90.0) % 360.0 array_rotation_rad = np.radians(array_rotation) + if ( + params.pv_array.torque_tube_separation > 0.0 + and params.pv_array.torque_tube_outer_radius > 0.0 + ): + self.modeling_torque_tube = True + else: + self.modeling_torque_tube = False + assert ( + params.pv_array.modules_per_span == 1 + ), "When not modeling torque tube, modules_per_span must be 1." + # The centroid of each panel in the x-direction (these should start at x=0) x_centers = np.linspace( 0.0, @@ -645,257 +1339,727 @@ def Rz(theta): domain_tag_list = [] # domain_tag_list.append(domain_tag) + # panel_tag_list includes all structure, panel+connector panel_tag_list = [] panel_ct = 0 - panel_id_y = -1 + # only include panels + structure_panel_only_list = [] + # only include connectors + structure_connector_only_list = [] - prev_surf_tag = [] - for k, yy in enumerate(y_centers): - panel_id_x = -1 - panel_id_y += 1 - for j, xx in enumerate(x_centers): - panel_id_x += 1 + module_distances = np.linspace( + -half_span, half_span, params.pv_array.modules_per_span + 1 + ) + + transformed_com = {} # all panels + + # panel_id_y = -1 + + # prev_surf_tag = [] + for panel_id_y, yy in enumerate(y_centers): + # panel_id_x = -1 + # panel_id_y += 1 + for panel_id_x, xx in enumerate(x_centers): + # panel_id_x += 1 # Create an 0-tracking-degree panel centered at (x, y, z) = (0, 0, 0) - panel_id = self.gmsh_model.occ.addBox( - -half_chord, - -half_span, - -half_thickness, - params.pv_array.panel_chord, - params.pv_array.panel_span, - params.pv_array.panel_thickness, - ) - panel_tag = (self.ndim, panel_id) - panel_tag_list.append(panel_tag) + this_panel_tag_list = [] # for each row + this_panel_transformed_com = {} # for each row - numpy_pt_list = [] - embedded_lines_tag_list = [] + numpy_pt_list = [] # for each row + embedded_lines_tag_list = [] # for each row - # Add a bisecting line to the bottom of the panel in the spanwise direction - pt_1 = self.gmsh_model.occ.addPoint(0, -half_span, -half_thickness) - pt_2 = self.gmsh_model.occ.addPoint(0, half_span, -half_thickness) + for module_id in range(params.pv_array.modules_per_span): - numpy_pt_list.append( - [0, -half_span, -half_thickness, 0, half_span, -half_thickness] - ) + module_span = ( + module_distances[module_id + 1] - module_distances[module_id] + ) - torque_tube_id = self.gmsh_model.occ.addLine(pt_1, pt_2) - torque_tube_tag = (1, torque_tube_id) - embedded_lines_tag_list.append(torque_tube_tag) - - # Add lines in the streamwise direction to mimic sections of panel held rigid by motor - if params.pv_array.span_fixation_pts is not None: - if not isinstance(params.pv_array.span_fixation_pts, list): - num_fixation_pts = int( - np.floor( - params.pv_array.panel_span - / params.pv_array.span_fixation_pts - ) + if ( + self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + ): + # Create an 0-tracking-degree panel centered at (x, y, z) = (0, 0, 0) + this_module = self.gmsh_model.occ.addBox( + -half_chord, + module_distances[module_id], + params.pv_array.torque_tube_separation - half_thickness, + params.pv_array.panel_chord, + module_span, + params.pv_array.panel_thickness, ) - fixation_pts_list = [] - - for k in range(1, num_fixation_pts + 1): - next_pt = k * params.pv_array.span_fixation_pts + panel_tag_list.append((self.ndim, this_module)) # all panels + this_panel_tag_list.append( + (self.ndim, this_module) + ) # this panel array/row + + structure_panel_only_list.append((self.ndim, this_module)) + + this_standoff = self.gmsh_model.occ.addBox( + -params.pv_array.block_chord_div_by_panel_chord + * half_chord, + module_distances[module_id], + -half_thickness, + 2.0 + * params.pv_array.block_chord_div_by_panel_chord + * half_chord, + params.pv_array.block_chord_div_by_panel_chord * half_chord, + params.pv_array.torque_tube_separation, + ) + panel_tag_list.append((self.ndim, this_standoff)) + this_panel_tag_list.append((self.ndim, this_standoff)) - eps = 1e-9 + structure_connector_only_list.append((self.ndim, this_standoff)) - if ( - next_pt > eps - and next_pt < params.pv_array.panel_span - eps - ): - fixation_pts_list.append(next_pt) + # Add a bisecting line to the bottom of the connector in the spanwise direction + pt_1 = self.gmsh_model.occ.addPoint( + 0, module_distances[module_id], -half_thickness + ) + pt_2 = self.gmsh_model.occ.addPoint( + 0, + module_distances[module_id] + + params.pv_array.block_chord_div_by_panel_chord + * half_chord, + -half_thickness, + ) + numpy_pt_list.append( + [ + 0, + module_distances[module_id], + -half_thickness, + 0, + module_distances[module_id] + + params.pv_array.block_chord_div_by_panel_chord + * half_chord, + -half_thickness, + ] + ) # for this row + torque_tube_id = self.gmsh_model.occ.addLine(pt_1, pt_2) + torque_tube_tag = (1, torque_tube_id) + embedded_lines_tag_list.append( + torque_tube_tag + ) # for this panel row else: - fixation_pts_list = params.pv_array.span_fixation_pts + this_module = self.gmsh_model.occ.addBox( + -half_chord, + module_distances[module_id], + -half_thickness, + params.pv_array.panel_chord, + module_span, + params.pv_array.panel_thickness, + ) + panel_tag_list.append((self.ndim, this_module)) + this_panel_tag_list.append((self.ndim, this_module)) + + structure_panel_only_list.append((self.ndim, this_module)) - for fp in fixation_pts_list: pt_1 = self.gmsh_model.occ.addPoint( - -half_chord, -half_span + fp, -half_thickness + 0, module_distances[module_id], -half_thickness ) pt_2 = self.gmsh_model.occ.addPoint( - half_chord, -half_span + fp, -half_thickness + 0, module_distances[module_id + 1], -half_thickness ) - - # FIXME: don't add the fixation points into the numpy tagging for now numpy_pt_list.append( [ - -half_chord, - -half_span + fp, + 0, + module_distances[module_id], -half_thickness, - half_chord, - -half_span + fp, + 0, + module_distances[module_id + 1], -half_thickness, ] ) + torque_tube_id = self.gmsh_model.occ.addLine(pt_1, pt_2) + torque_tube_tag = (1, torque_tube_id) + embedded_lines_tag_list.append( + torque_tube_tag + ) # for this panel row + + # Add lines in the streamwise direction to mimic sections of panel held rigid by motor + if params.pv_array.span_fixation_pts is not None: + if not isinstance(params.pv_array.span_fixation_pts, list): + num_fixation_pts = int( + np.floor( + params.pv_array.panel_span + / params.pv_array.span_fixation_pts + ) + ) + + fixation_pts_list = [] + + for k in range(1, num_fixation_pts + 1): + next_pt = k * params.pv_array.span_fixation_pts + + eps = 1e-9 + + if ( + next_pt > eps + and next_pt < params.pv_array.panel_span - eps + ): + fixation_pts_list.append(next_pt) + + else: + fixation_pts_list = params.pv_array.span_fixation_pts + + for fp in fixation_pts_list: + pt_1 = self.gmsh_model.occ.addPoint( + -half_chord, -half_span + fp, -half_thickness + ) + pt_2 = self.gmsh_model.occ.addPoint( + half_chord, -half_span + fp, -half_thickness + ) + + # FIXME: don't add the fixation points into the numpy tagging for now + numpy_pt_list.append( + [ + -half_chord, + -half_span + fp, + -half_thickness, + half_chord, + -half_span + fp, + -half_thickness, + ] + ) + + fixed_pt_id = self.gmsh_model.occ.addLine(pt_1, pt_2) + fixed_pt_tag = (1, fixed_pt_id) + + embedded_lines_tag_list.append(fixed_pt_tag) + + if ( + self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + ): # add the last connector + last_standoff = self.gmsh_model.occ.addBox( + -params.pv_array.block_chord_div_by_panel_chord * half_chord, + module_distances[module_id + 1] + - params.pv_array.block_chord_div_by_panel_chord * half_chord, + -half_thickness, + 2.0 + * params.pv_array.block_chord_div_by_panel_chord + * half_chord, + params.pv_array.block_chord_div_by_panel_chord * half_chord, + params.pv_array.torque_tube_separation, + ) + this_panel_tag_list.append((self.ndim, last_standoff)) + panel_tag_list.append((self.ndim, last_standoff)) - fixed_pt_id = self.gmsh_model.occ.addLine(pt_1, pt_2) - fixed_pt_tag = (1, fixed_pt_id) + structure_connector_only_list.append((self.ndim, last_standoff)) - embedded_lines_tag_list.append(fixed_pt_tag) + pt_1 = self.gmsh_model.occ.addPoint( + 0, + module_distances[module_id + 1] + - params.pv_array.block_chord_div_by_panel_chord * half_chord, + -half_thickness, + ) + pt_2 = self.gmsh_model.occ.addPoint( + 0, module_distances[module_id + 1], -half_thickness + ) + numpy_pt_list.append( + [ + 0, + module_distances[module_id + 1] + - params.pv_array.block_chord_div_by_panel_chord + * half_chord, + -half_thickness, + 0, + module_distances[module_id + 1], + -half_thickness, + ] + ) + torque_tube_id = self.gmsh_model.occ.addLine(pt_1, pt_2) + torque_tube_tag = (1, torque_tube_id) + embedded_lines_tag_list.append(torque_tube_tag) # Store the result of fragmentation, it holds all the small surfaces we need to tag panel_frags = self.gmsh_model.occ.fragment( - [panel_tag], embedded_lines_tag_list - ) - - # Translate the panel by (x_center, y_center, elev) - self.gmsh_model.occ.translate( - [panel_tag], - xx, - yy, - params.pv_array.elevation, - ) - - # extract just the first entry, and remove the 3d entry in position 0 - panel_surfs = panel_frags[0] - panel_surfs.pop(0) - panel_surfs = [k[1] for k in panel_surfs] - - # TODO: USE THESE UNAMBIGUOUS NAMES IN A FUTURE REFACTOR - # self._add_to_domain_markers(f"x_min_{panel_ct:.0f}", [panel_surfs[0]], "facet") - # self._add_to_domain_markers(f"x_max_{panel_ct:.0f}", [panel_surfs[1]], "facet") - # self._add_to_domain_markers(f"y_min_{panel_ct:.0f}", [panel_surfs[2]], "facet") - # self._add_to_domain_markers(f"y_max_{panel_ct:.0f}", [panel_surfs[3]], "facet") - # self._add_to_domain_markers(f"z_min_{panel_ct:.0f}", panel_surfs[4:-1], "facet") - # self._add_to_domain_markers(f"z_max_{panel_ct:.0f}", [panel_surfs[-1]], "facet") - - # self._add_to_domain_markers(f"front_{panel_ct:.0f}", [panel_surfs[0]], "facet") - # self._add_to_domain_markers(f"back_{panel_ct:.0f}", [panel_surfs[1]], "facet") - # self._add_to_domain_markers(f"left_{panel_ct:.0f}", [panel_surfs[2]], "facet") - # self._add_to_domain_markers(f"right_{panel_ct:.0f}", [panel_surfs[3]], "facet") - # self._add_to_domain_markers(f"bottom_{panel_ct:.0f}", panel_surfs[4:-1], "facet") - # self._add_to_domain_markers(f"top_{panel_ct:.0f}", [panel_surfs[-1]], "facet") - - # self._add_to_domain_markers( - # f"front_{panel_ct:.0f}", [panel_surfs[1]], "facet" - # ) # should be bottom - # self._add_to_domain_markers( - # f"back_{panel_ct:.0f}", [panel_surfs[2]], "facet" - # ) - - # self._add_to_domain_markers( - # f"left_{panel_ct:.0f}", [panel_surfs[3]], "facet" - # ) # should be left - # self._add_to_domain_markers( - # f"right_{panel_ct:.0f}", [panel_surfs[4]], "facet" - # ) # correct 4 - - # self._add_to_domain_markers( - # f"bottom_{panel_ct:.0f}", panel_surfs[5:-1], "facet" - # ) - # self._add_to_domain_markers( - # f"top_{panel_ct:.0f}", [panel_surfs[-1]], "facet" - # ) # should be front - - top_coord = ( - params.domain.z_min - + (params.pv_array.elevation - params.domain.z_min) - + params.pv_array.panel_thickness / 2 - ) - print("top", top_coord) - bottom_coord = ( - params.domain.z_min - + (params.pv_array.elevation - params.domain.z_min) - - params.pv_array.panel_thickness / 2 + this_panel_tag_list, embedded_lines_tag_list ) - print("bottom", bottom_coord) - left_coord = -params.pv_array.panel_chord / 2 + panel_id_x * ( - params.pv_array.stream_spacing - ) - print("left", left_coord) - right_coord = +params.pv_array.panel_chord / 2 + panel_id_x * ( - params.pv_array.stream_spacing - ) - print("right", right_coord) - - front_coord = -params.pv_array.panel_span / 2 + yy - print("front", front_coord) - back_coord = +params.pv_array.panel_span / 2 + yy - print("back", back_coord) - surf_tag_list_total = self.gmsh_model.occ.getEntities(self.ndim - 1) + for panel_tag in this_panel_tag_list: # 3d cell domain + self.gmsh_model.occ.synchronize() - surf_tag_list = [ - vector - for vector in surf_tag_list_total - if vector not in prev_surf_tag - ] + # Get the list of 2D surfaces (surfaces) that make up this panel + surf_tags_for_this_panel = self.gmsh_model.getBoundary( + [panel_tag], oriented=False + ) - prev_surf_tag = surf_tag_list_total - for surf_tag in surf_tag_list: - surf_id = surf_tag[1] - com = self.gmsh_model.occ.getCenterOfMass(self.ndim - 1, surf_id) - print(com) - # sturctures tagging - if np.isclose(com[2], bottom_coord): - self._add_to_domain_markers( - f"bottom_{panel_ct:.0f}", [surf_id], "facet" - ) - print("bottom found") - # self._add_to_domain_markers("x_min", [surf_id], "facet") + vol_com = self.gmsh_model.occ.getCenterOfMass( + self.ndim, panel_tag[1] + ) + if np.isclose( + vol_com[2], + params.pv_array.torque_tube_separation, + ): + target_key = f"modules" + elif np.isclose( + vol_com[2], + params.pv_array.torque_tube_separation / 2.0 - half_thickness, + ): + target_key = f"connectors" + + if target_key is not None: + if target_key in this_panel_transformed_com: + this_panel_transformed_com[target_key].append(vol_com) + else: + this_panel_transformed_com[target_key] = [vol_com] + + for surf_tag in surf_tags_for_this_panel: + surf_dim = surf_tag[0] + surf_id = surf_tag[1] + com = self.gmsh_model.occ.getCenterOfMass(surf_dim, surf_id) + + target_key = None + + surface_located_or_not = False + + # sturctures tagging + if np.isclose(com[0], -half_chord): + target_key = f"panel_left_{panel_ct:.0f}" + surface_located_or_not = True + + if np.isclose(com[0], half_chord): + target_key = f"panel_right_{panel_ct:.0f}" + surface_located_or_not = True + + if ( + np.isclose(com[1], module_distances[0]) + and np.isclose(com[0], 0) + and np.isclose( + com[2], + params.pv_array.torque_tube_separation, + ) + ): + target_key = f"panel_front_{panel_ct:.0f}" + surface_located_or_not = True + + if ( + np.isclose( + com[1], + module_distances[params.pv_array.modules_per_span], + ) + and np.isclose(com[0], 0) + and np.isclose( + com[2], + params.pv_array.torque_tube_separation, + ) + ): + target_key = f"panel_back_{panel_ct:.0f}" + surface_located_or_not = True + + if (not self.modeling_torque_tube) or ( + params.general.geometry_module != "panels3d" + ): + if np.isclose( + com[2], + params.pv_array.torque_tube_separation - half_thickness, + ): + target_key = f"panel_bottom_{panel_ct:.0f}_0" + surface_located_or_not = True - elif np.allclose(com[2], top_coord): - self._add_to_domain_markers( - f"top_{panel_ct:.0f}", [surf_id], "facet" - ) - print("top found") - # self._add_to_domain_markers("x_max", [surf_id], "facet") + if np.isclose( + com[2], + half_thickness + params.pv_array.torque_tube_separation, + ): + target_key = f"panel_top_{panel_ct:.0f}_0" + surface_located_or_not = True + + else: + for module_id in range(params.pv_array.modules_per_span): + if ( + com[1] + >= module_distances[module_id] + + module_span / 2.0 + - params.pv_array.block_chord_div_by_panel_chord + * half_chord + and com[1] + <= module_distances[module_id] + + module_span / 2.0 + + params.pv_array.block_chord_div_by_panel_chord + * half_chord + and np.isclose(com[0], 0) + and np.isclose( + com[2], + params.pv_array.torque_tube_separation + - half_thickness, + ) + ): + target_key = ( + f"panel_bottom_{panel_ct:.0f}_{module_id:.0f}" + ) + surface_located_or_not = True + break + + if ( + com[1] + >= module_distances[module_id] + + module_span / 2.0 + - params.pv_array.block_chord_div_by_panel_chord + * half_chord + and com[1] + <= module_distances[module_id] + + module_span / 2.0 + + params.pv_array.block_chord_div_by_panel_chord + * half_chord + and np.isclose(com[0], 0) + and np.isclose( + com[2], + params.pv_array.torque_tube_separation + + half_thickness, + ) + ): + target_key = ( + f"panel_top_{panel_ct:.0f}_{module_id:.0f}" + ) + surface_located_or_not = True + break + # if ( + # np.isclose( + # com[2], + # params.pv_array.torque_tube_separation + # + params.pv_array.panel_thickness + # ) + + # ): + # target_key = f"panel_top_{panel_ct:.0f}" + # surface_located_or_not = True + + # mark the last block in each row + if ( + np.isclose( + com[0], + params.pv_array.block_chord_div_by_panel_chord + * half_chord, + ) + and self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + and np.isclose( + com[1], + module_distances[params.pv_array.modules_per_span] + - params.pv_array.block_chord_div_by_panel_chord + * half_chord + / 2.0, + ) + ): + target_key = f"block_right_{panel_ct:.0f}_{params.pv_array.modules_per_span:.0f}" + surface_located_or_not = True - elif np.allclose(com[0], left_coord): - self._add_to_domain_markers( - f"left_{panel_ct:.0f}", [surf_id], "facet" - ) - print("left found") - # self._add_to_domain_markers("y_min", [surf_id], "facet") + if ( + np.isclose( + com[0], + -params.pv_array.block_chord_div_by_panel_chord + * half_chord, + ) + and self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + and np.isclose( + com[1], + module_distances[params.pv_array.modules_per_span] + - params.pv_array.block_chord_div_by_panel_chord + * half_chord + / 2.0, + ) + ): + target_key = f"block_left_{panel_ct:.0f}_{params.pv_array.modules_per_span:.0f}" + surface_located_or_not = True - elif np.allclose(com[0], right_coord): - self._add_to_domain_markers( - f"right_{panel_ct:.0f}", [surf_id], "facet" - ) - print("right found") + if ( + np.isclose( + com[2], + params.pv_array.torque_tube_separation / 2.0 + - half_thickness, + ) + and self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + and np.isclose( + com[1], + module_distances[params.pv_array.modules_per_span] + - params.pv_array.block_chord_div_by_panel_chord + * half_chord, + ) + ): + target_key = f"block_front_{panel_ct:.0f}_{params.pv_array.modules_per_span:.0f}" + surface_located_or_not = True - elif np.allclose(com[1], front_coord): - self._add_to_domain_markers( - f"front_{panel_ct:.0f}", [surf_id], "facet" - ) - print("front found") - # self._add_to_domain_markers("y_min", [surf_id], "facet") + if ( + np.isclose( + com[2], + params.pv_array.torque_tube_separation / 2.0 + - half_thickness, + ) + and self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + and np.isclose( + com[1], + module_distances[params.pv_array.modules_per_span], + ) + ): + target_key = f"block_back_{panel_ct:.0f}_{params.pv_array.modules_per_span:.0f}" + surface_located_or_not = True - elif np.allclose(com[1], back_coord): - self._add_to_domain_markers( - f"back_{panel_ct:.0f}", [surf_id], "facet" - ) - print("back found") - # self._add_to_domain_markers("y_max", [surf_id], "facet") + if ( + np.isclose(com[2], -half_thickness) + and self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + and np.isclose( + com[1], + module_distances[params.pv_array.modules_per_span] + - params.pv_array.block_chord_div_by_panel_chord + * half_chord + / 2.0, + ) + ): + target_key = f"block_bottom_{panel_ct:.0f}_{params.pv_array.modules_per_span:.0f}" + surface_located_or_not = True + if ( + np.isclose( + com[2], + params.pv_array.torque_tube_separation + - half_thickness, + ) + and self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + and np.isclose( + com[1], + module_distances[params.pv_array.modules_per_span] + - params.pv_array.block_chord_div_by_panel_chord + * half_chord + / 2.0, + ) + ): + target_key = f"interior_surface_{panel_ct:.0f}" # block/panel interface + surface_located_or_not = True + + # if not the most left/right panel boundary, it is panel/panel interface + if not surface_located_or_not: + for module_id in range( + params.pv_array.modules_per_span + ): + + # sturctures tagging + + if ( + np.isclose(com[1], module_distances[module_id]) + and np.isclose(com[0], 0) + and np.isclose( + com[2], + params.pv_array.torque_tube_separation, + ) + ): + target_key = f"interior_surface_{panel_ct:.0f}" + surface_located_or_not = True + break + + if ( + np.isclose( + com[1], module_distances[module_id + 1] + ) + and np.isclose(com[0], 0) + and np.isclose( + com[2], + params.pv_array.torque_tube_separation, + ) + ): + target_key = f"interior_surface_{panel_ct:.0f}" + surface_located_or_not = True + break + + # if ( + # np.isclose(com[2], params.pv_array.torque_tube_separation) + # and np.isclose(com[0], 0.0) and com[1] >= module_distances[module_id]+module_span/2.0-params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0 + # and com[1] <= module_distances[module_id]+module_span/2.0+params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0 + # ): + # target_key = f"panel_bottom_{panel_ct:.0f}" + # surface_located_or_not = True + # break + + if ( + np.isclose( + com[2], + params.pv_array.torque_tube_separation + - half_thickness, + ) + and np.isclose(com[0], 0.0) + and np.isclose( + com[1], + module_distances[module_id] + + params.pv_array.block_chord_div_by_panel_chord + * half_chord + / 2.0, + ) + and self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + ): + target_key = f"interior_surface_{panel_ct:.0f}" + surface_located_or_not = True + break + + if ( + np.isclose( + com[0], + params.pv_array.block_chord_div_by_panel_chord + * half_chord, + ) + and self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + and np.isclose( + com[1], + module_distances[module_id] + + params.pv_array.block_chord_div_by_panel_chord + * half_chord + / 2.0, + ) + ): + target_key = f"block_right_{panel_ct:.0f}_{module_id:.0f}" + surface_located_or_not = True + break + + if ( + np.isclose( + com[0], + -params.pv_array.block_chord_div_by_panel_chord + * half_chord, + ) + and self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + and np.isclose( + com[1], + module_distances[module_id] + + params.pv_array.block_chord_div_by_panel_chord + * half_chord + / 2.0, + ) + ): + target_key = ( + f"block_left_{panel_ct:.0f}_{module_id:.0f}" + ) + surface_located_or_not = True + break + + if ( + np.isclose( + com[2], + params.pv_array.torque_tube_separation / 2.0 + - half_thickness, + ) + and self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + and np.isclose( + com[1], module_distances[module_id] + ) + ): + target_key = f"block_front_{panel_ct:.0f}_{module_id:.0f}" + surface_located_or_not = True + break + + if ( + np.isclose( + com[2], + params.pv_array.torque_tube_separation / 2.0 + - half_thickness, + ) + and self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + and np.isclose( + com[1], + module_distances[module_id] + + params.pv_array.block_chord_div_by_panel_chord + * half_chord, + ) + ): + target_key = ( + f"block_back_{panel_ct:.0f}_{module_id:.0f}" + ) + surface_located_or_not = True + break + + if ( + np.isclose(com[2], -half_thickness) + and self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + and np.isclose( + com[1], + module_distances[module_id] + + params.pv_array.block_chord_div_by_panel_chord + * half_chord + / 2.0, + ) + ): + target_key = f"block_bottom_{panel_ct:.0f}_{module_id:.0f}" + surface_located_or_not = True + break + + if not surface_located_or_not: + target_key = f"trash_{panel_ct:.0f}" + print("facet in trash") + + if target_key is not None: + if target_key in this_panel_transformed_com: + this_panel_transformed_com[target_key].append(com) + else: + this_panel_transformed_com[target_key] = [com] + + # print(this_panel_transformed_com[f"trash_{panel_ct:.0f}"]) + + for ( + key, + val, + ) in ( + this_panel_transformed_com.items() + ): # all facets for each panel row. + for row_num, com in enumerate(val): + com_array = np.array(com) + + com_array = np.dot(com_array, Ry(tracker_angle_rad).T) + + com_array[0] += xx + com_array[1] += yy + com_array[2] += params.pv_array.elevation + + com_array[0] -= x_center_of_mass + com_array[1] -= y_center_of_mass + + com_array = np.dot(com_array, Rz(array_rotation_rad).T) + + com_array[0] += x_center_of_mass + com_array[1] += y_center_of_mass + + if key in transformed_com: + transformed_com[key].append(com_array) + else: + transformed_com[key] = [com_array] # including all rows + + # actually this is panel array count panel_ct += 1 - # Rotate the panel by its tracking angle along the y-axis (currently centered at (0, 0, 0)) + # Rotate the panel by its tracking angle along the y-axis + # (currently centered at (0.0, 0.0, 0.0)) self.gmsh_model.occ.rotate( - [panel_tag], 0, 0, 0, 0, 1, 0, tracker_angle_rad + this_panel_tag_list, + 0.0, + 0.0, + 0.0, + 0, + 1, + 0, + tracker_angle_rad, ) - numpy_pt_panel_array = np.array(numpy_pt_list) - numpy_pt_panel_array = np.reshape(numpy_pt_panel_array, (-1, self.ndim)) - - numpy_pt_panel_array = np.dot( - numpy_pt_panel_array, Ry(tracker_angle_rad).T + # Translate the panel by (x_center, y_center, elev) + self.gmsh_model.occ.translate( + this_panel_tag_list, + xx, + yy, + params.pv_array.elevation, ) - # if not hasattr(self, "numpy_pt_array"): - # numpy_pt_array = np.array(numpy_pt_list) - # else: - # numpy_pt_array = np.vcat(numpy_pt_array, np.array(numpy_pt_list)) - - numpy_pt_panel_array[:, 0] += xx - numpy_pt_panel_array[:, 1] += yy - numpy_pt_panel_array[:, 2] += params.pv_array.elevation - # Rotate the panel about the center of the full array as a proxy for changing wind direction (x_center, y_center, 0) self.gmsh_model.occ.rotate( - [panel_tag], + this_panel_tag_list, x_center_of_mass, y_center_of_mass, 0, @@ -905,6 +2069,22 @@ def Rz(theta): array_rotation_rad, ) + # Now, apply the same transformations to the numpy representation + # Rotate the panel by its tracking angle along the y-axis + # (currently centered at (0.0, 0.0, 0.0)) + numpy_pt_panel_array = np.array(numpy_pt_list) + numpy_pt_panel_array = np.reshape(numpy_pt_panel_array, (-1, self.ndim)) + + numpy_pt_panel_array = np.dot( + numpy_pt_panel_array, Ry(tracker_angle_rad).T + ) + + # Translate the panel by (x_center, y_center, elev) + numpy_pt_panel_array[:, 0] += xx + numpy_pt_panel_array[:, 1] += yy + numpy_pt_panel_array[:, 2] += params.pv_array.elevation + + # Rotate the panel about the center of the full array as a proxy for changing wind direction (x_center, y_center, 0) numpy_pt_panel_array[:, 0] -= x_center_of_mass numpy_pt_panel_array[:, 1] -= y_center_of_mass @@ -922,85 +2102,84 @@ def Rz(theta): else: self.numpy_pt_total_array = np.copy(numpy_pt_panel_array) - # Check that this panel still exists in the confines of the domain - bbox = self.gmsh_model.occ.get_bounding_box(panel_tag[0], panel_tag[1]) - - if bbox[0] < params.domain.x_min: - raise ValueError( - f"Panel with location (x, y) = ({xx}, {yy}) extends past x_min wall." - ) - if bbox[1] < params.domain.y_min: - raise ValueError( - f"Panel with location (x, y) = ({xx}, {yy}) extends past y_min wall." - ) - if bbox[3] > params.domain.x_max: - raise ValueError( - f"Panel with location (x, y) = ({xx}, {yy}) extends past x_max wall." - ) - if bbox[4] > params.domain.y_max: - raise ValueError( - f"Panel with location (x, y) = ({xx}, {yy}) extends past y_max wall." - ) - # Fragment all panels from the overall domain self.gmsh_model.occ.fragment(domain_tag_list, panel_tag_list) self.gmsh_model.occ.synchronize() + gmsh.write("pnales.brep") + # exit() self.numpy_pt_total_array = np.reshape( self.numpy_pt_total_array, (-1, int(2 * self.ndim)) ) - # import matplotlib.pyplot as plt - # for k in self.numpy_pt_total_array: - # plt.plot([k[0], k[3]], [k[1], k[4]]) - # plt.show() - - # # Surfaces are the entities with dimension equal to the mesh dimension -1 - # surf_tag_list = self.gmsh_model.occ.getEntities(self.ndim-1) - - # for surf_tag in surf_tag_list: - # surf_id = surf_tag[1] - # com = self.gmsh_model.occ.getCenterOfMass(self.ndim-1, surf_id) - - # #sturctures tagging - # if np.isclose(com[0], params.domain.x_min): - # self._add_to_domain_markers("x_min", [surf_id], "facet") + # print(transformed_com) + # exit() - # elif np.allclose(com[0], params.domain.x_max): - # self._add_to_domain_markers("x_max", [surf_id], "facet") + # for all panels + # Loop over all the finalized surfaces after fragmentation and tag everything + all_surf_tag_list = self.gmsh_model.occ.getEntities(self.ndim - 1) - # elif np.allclose(com[1], params.domain.y_min): - # self._add_to_domain_markers("y_min", [surf_id], "facet") + for surf_tag in all_surf_tag_list: + surf_id = surf_tag[1] + com = self.gmsh_model.occ.getCenterOfMass(self.ndim - 1, surf_id) - # elif np.allclose(com[1], params.domain.y_max): - # self._add_to_domain_markers("y_max", [surf_id], "facet") + located_this_surface = False + + for key, val in transformed_com.items(): + for target_com in val: + # print(target_com) + if np.allclose(np.array(com), target_com): + located_this_surface = True + if "trash" not in key: + # print(key) + self._add_to_domain_markers(key, [surf_id], "facet") + # if "interior_surface" not in key: + # self._add_to_domain_markers("structure_fluid_interface", [surf_id], "facet") + + if not located_this_surface: + print( + f"Warning: Surface {surf_tag} has not been added to domain markers" + ) - # elif np.allclose(com[2], params.domain.z_min): - # self._add_to_domain_markers("z_min", [surf_id], "facet") + # Since this is not one of the exterior walls, we should check if it extends + # past the boundaries x_min, x_max, ... + this_surf_bbox = self.gmsh_model.occ.get_bounding_box( + self.ndim - 1, surf_id + ) - # elif np.allclose(com[2], params.domain.z_max): - # self._add_to_domain_markers("z_max", [surf_id], "facet") + # Test that the rotated point still exists in the box domain + if this_surf_bbox[0] < params.domain.x_min: + raise ValueError(f"A panel extends past the x_min wall.") + if this_surf_bbox[0] > params.domain.x_max: + raise ValueError(f"A panel extends past the x_max wall.") + if this_surf_bbox[1] < params.domain.y_min: + raise ValueError(f"A panel extends past the y_min wall.") + if this_surf_bbox[1] > params.domain.y_max: + raise ValueError(f"A panel extends past the y_max wall.") + if this_surf_bbox[2] < 0.0: + raise ValueError( + f"A panel extends past the z_min wall (ground level = 0.0)." + ) + if this_surf_bbox[2] > params.domain.z_max: + raise ValueError(f"A panel extends past the z_max wall.") - # Volumes are the entities with dimension equal to the mesh dimension - vol_tag_list = self.gmsh_model.occ.getEntities(self.ndim) - structure_vol_list = [] - # fluid_vol_list = [] + # mark the panel and connector volumes + all_vol_tag_list = self.gmsh_model.occ.getEntities(self.ndim) - for vol_tag in vol_tag_list: + for vol_tag in all_vol_tag_list: vol_id = vol_tag[1] - structure_vol_list.append(vol_id) - # vol_id = vol_tag[1] - - # if vol_id <= params.pv_array.stream_rows * params.pv_array.span_rows: - # # Solid Cell + com = self.gmsh_model.occ.getCenterOfMass(self.ndim, vol_id) - # else: - # # Fluid Cell - # fluid_vol_list.append(vol_id) + located_this_volume = False - self._add_to_domain_markers("structure", structure_vol_list, "cell") - # self._add_to_domain_markers("fluid", fluid_vol_list, "cell") + for key, val in transformed_com.items(): + for target_com in val: + # print(target_com) + if np.allclose(np.array(com), target_com): + located_this_volume = True + if "trash" not in key: + self._add_to_domain_markers(key, [vol_id], "cell") # Record all the data collected to domain_markers as physics groups with physical names # because we are creating domain_markers _within_ the build method, we don't need to @@ -1144,25 +2323,47 @@ def set_length_scales(self, params, domain_markers): internal_surface_tags = [] for panel_id in range(params.pv_array.stream_rows * params.pv_array.span_rows): - internal_surface_tags.append( - domain_markers[f"bottom_{panel_id}"]["gmsh_tags"][0] - ) - internal_surface_tags.append( - domain_markers[f"top_{panel_id}"]["gmsh_tags"][0] + + internal_surface_tags.extend( + domain_markers[f"panel_left_{panel_id}"]["gmsh_tags"] ) - internal_surface_tags.append( - domain_markers[f"left_{panel_id}"]["gmsh_tags"][0] + internal_surface_tags.extend( + domain_markers[f"panel_right_{panel_id}"]["gmsh_tags"] ) - internal_surface_tags.append( - domain_markers[f"right_{panel_id}"]["gmsh_tags"][0] + internal_surface_tags.extend( + domain_markers[f"panel_front_{panel_id}"]["gmsh_tags"] ) - internal_surface_tags.append( - domain_markers[f"front_{panel_id}"]["gmsh_tags"][0] - ) - internal_surface_tags.append( - domain_markers[f"back_{panel_id}"]["gmsh_tags"][0] + internal_surface_tags.extend( + domain_markers[f"panel_back_{panel_id}"]["gmsh_tags"] ) + if ( + self.modeling_torque_tube + and params.general.geometry_module == "panels3d" + ): + for module_id in range(params.pv_array.modules_per_span): + internal_surface_tags.extend( + domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"][ + "gmsh_tags" + ] + ) + internal_surface_tags.extend( + domain_markers[f"panel_top_{panel_id:.0f}_{module_id:.0f}"][ + "gmsh_tags" + ] + ) + else: # delete after testing + internal_surface_tags.extend( + domain_markers[f"panel_bottom_{panel_id}_0"]["gmsh_tags"] + ) + internal_surface_tags.extend( + domain_markers[f"panel_top_{panel_id}_0"]["gmsh_tags"] + ) + + # print(internal_surface_tags) + # print(len(internal_surface_tags)) + # exit() + min_dist = [] # Define a distance field from the immersed panels diff --git a/pvade/structure/ElasticityAnalysis.py b/pvade/structure/ElasticityAnalysis.py index 4a84d988..fe320d22 100644 --- a/pvade/structure/ElasticityAnalysis.py +++ b/pvade/structure/ElasticityAnalysis.py @@ -87,6 +87,381 @@ def build_boundary_conditions(self, domain, params): self.bc = build_structure_boundary_conditions(domain, params, self.V) + def calculate_K_for_Robin_BC(self, domain, flow, params): + + # shear modulus of tube + Gt = params.structure.elasticity_modulus_tube / ( + 2 * (1 + params.structure.poissons_ratio_tube) + ) + + Ipt = ( + np.pi + / 2 + * ( + params.pv_array.torque_tube_outer_radius**4 + - params.pv_array.torque_tube_inner_radius**4 + ) + ) + + Gp = params.structure.elasticity_modulus / ( + 2 * (1 + params.structure.poissons_ratio) + ) + + Ipp = 1 / 3.0 * params.pv_array.panel_thickness * params.pv_array.panel_chord**3 + + # total number of pv rows + total_num_panels = params.pv_array.stream_rows * params.pv_array.span_rows + + for panel_id in range(total_num_panels): + # construct the matrix, size: modules_per_span+1 by modules_per_span+1 + theo_matrix = np.zeros( + ( + params.pv_array.modules_per_span + 1, + params.pv_array.modules_per_span + 1, + ) + ) + theo_vector = np.zeros(params.pv_array.modules_per_span + 1) + + connector_locations = np.linspace( + 0, params.pv_array.panel_span, params.pv_array.modules_per_span + 1 + ) + + # ?? double integral is not available until the second time step. Should make it be zero for the + # first time step then it will be updated once the first step is finished? but the assemble is out of + # the time loop. Should we solve fluid then structure? + + total_torque_on_this_panel_name = f"total_torque_panel_{panel_id:.0f}" + total_torque_on_this_panel = getattr( + flow, total_torque_on_this_panel_name + ).value + + theo_matrix[params.pv_array.modules_per_span, :] = 1 / Gp / Ipp + theo_vector[params.pv_array.modules_per_span] = ( + -total_torque_on_this_panel / Gp / Ipp + ) + + # contribution from panel: + + # contribution from C0 to matrix: + + theo_matrix[ + : params.pv_array.modules_per_span, : params.pv_array.fixed_location + ] += connector_locations[params.pv_array.fixed_location] / (Gp * Ipp) + theo_matrix[ + : params.pv_array.modules_per_span, 1 : params.pv_array.fixed_location + ] -= connector_locations[1 : params.pv_array.fixed_location] / (Gp * Ipp) + + # contribution from C0 to vector: + T_double_integral_at_fixed_location_name = f"double_integral_total_torque_panel_{panel_id:.0f}_{params.pv_array.fixed_location:.0f}" + T_double_integral_at_fixed_location = getattr( + flow, T_double_integral_at_fixed_location_name + ).value + theo_vector[: params.pv_array.modules_per_span] += ( + -T_double_integral_at_fixed_location / Gp / Ipp + ) + + # contribution of panels to matrix + theo_matrix[0, 0] += -connector_locations[0] / (Gp * Ipp) + theo_matrix[1, 0] += -connector_locations[1] / (Gp * Ipp) + + for i in range(2, params.pv_array.fixed_location): + theo_matrix[i, :i] += -connector_locations[i] / (Gp * Ipp) + theo_matrix[i, 1:i] += connector_locations[1:i] / (Gp * Ipp) + + for i in range( + params.pv_array.fixed_location, params.pv_array.modules_per_span + ): + theo_matrix[i, : i + 1] += -connector_locations[i + 1] / (Gp * Ipp) + theo_matrix[i, 1 : i + 1] += connector_locations[1 : i + 1] / (Gp * Ipp) + + # contribution of panels to vector + T_double_integral_array = [] + for i in range(params.pv_array.modules_per_span + 1): + name = f"double_integral_total_torque_panel_{panel_id:.0f}_{i:.0f}" + T_double_integral_array.append(getattr(flow, name).value) + T_double_integral_array = np.array(T_double_integral_array) + theo_vector[: params.pv_array.modules_per_span] += ( + np.delete(T_double_integral_array, params.pv_array.fixed_location) + / Gp + / Ipp + ) + + # contribution from tube: + + # contribution from C0 to matrix: + + theo_matrix[ + : params.pv_array.modules_per_span, + params.pv_array.fixed_location : params.pv_array.modules_per_span + 1, + ] += -connector_locations[params.pv_array.fixed_location] / (Gt * Ipt) + theo_matrix[ + : params.pv_array.modules_per_span, 1 : params.pv_array.fixed_location + ] += -connector_locations[1 : params.pv_array.fixed_location] / (Gt * Ipt) + + # contribution from C0 to vector: + theo_vector[: params.pv_array.modules_per_span] += ( + total_torque_on_this_panel + * connector_locations[params.pv_array.fixed_location] + / Gt + / Ipt + ) + + # contribution of tube to matrix + theo_matrix[ + 0, 1 : params.pv_array.modules_per_span + 1 + ] += connector_locations[0] / (Gt * Ipt) + theo_matrix[ + 1, 1 : params.pv_array.modules_per_span + 1 + ] += connector_locations[1] / (Gt * Ipt) + + for i in range(2, params.pv_array.fixed_location): + theo_matrix[ + i, i : params.pv_array.modules_per_span + 1 + ] += connector_locations[i] / (Gt * Ipt) + theo_matrix[i, 1:i] += connector_locations[1:i] / (Gt * Ipt) + + for i in range( + params.pv_array.fixed_location, params.pv_array.modules_per_span + ): + theo_matrix[ + i, i + 1 : params.pv_array.modules_per_span + 1 + ] += connector_locations[i + 1] / (Gt * Ipt) + theo_matrix[i, 1 : i + 1] += connector_locations[1 : i + 1] / (Gt * Ipt) + + # contribution of tube to vector + theo_vector[: params.pv_array.fixed_location] += ( + -total_torque_on_this_panel + / Gt + / Ipt + * connector_locations[: params.pv_array.fixed_location] + ) + theo_vector[ + params.pv_array.fixed_location : params.pv_array.modules_per_span + ] += ( + -total_torque_on_this_panel + / Gt + / Ipt + * connector_locations[params.pv_array.fixed_location] + ) + + R_reaction_torque = np.dot( + np.linalg.inv(theo_matrix), theo_vector + ) # this is the torque applied to tube + + # rotation of tube at each connector location + tube_rotate_matrix = np.zeros( + (params.pv_array.modules_per_span, params.pv_array.modules_per_span + 1) + ) + tube_rotate_vector = np.zeros(params.pv_array.modules_per_span) + + # contribution from panel: + + # contribution from C0 to matrix: + + tube_rotate_matrix[ + : params.pv_array.modules_per_span, : params.pv_array.fixed_location + ] += connector_locations[params.pv_array.fixed_location] / (Gp * Ipp) + tube_rotate_matrix[ + : params.pv_array.modules_per_span, 1 : params.pv_array.fixed_location + ] -= connector_locations[1 : params.pv_array.fixed_location] / (Gp * Ipp) + + # contribution from C0 to vector: + tube_rotate_vector[: params.pv_array.modules_per_span] += ( + T_double_integral_array[params.pv_array.fixed_location] / Gp / Ipp + ) + + # contribution of panels to matrix + tube_rotate_matrix[0, 0] += -connector_locations[0] / (Gp * Ipp) + tube_rotate_matrix[1, 0] += -connector_locations[1] / (Gp * Ipp) + + for i in range(2, params.pv_array.fixed_location): + tube_rotate_matrix[i, :i] += -connector_locations[i] / (Gp * Ipp) + tube_rotate_matrix[i, 1:i] += connector_locations[1:i] / (Gp * Ipp) + + for i in range( + params.pv_array.fixed_location, params.pv_array.modules_per_span + ): + tube_rotate_matrix[i, : i + 1] += -connector_locations[i + 1] / ( + Gp * Ipp + ) + tube_rotate_matrix[i, 1 : i + 1] += connector_locations[1 : i + 1] / ( + Gp * Ipp + ) + + # contribution of panels to vector + tube_rotate_vector[: params.pv_array.modules_per_span] += ( + -np.delete(T_double_integral_array, params.pv_array.fixed_location) + / Gp + / Ipp + ) + + phi = np.dot(tube_rotate_matrix, R_reaction_torque) + tube_rotate_vector + + if isinstance(params.pv_array.tracker_angle, list): + if panel_id == 0: + assert ( + len(params.pv_array.tracker_angle) + == params.pv_array.stream_rows * params.pv_array.span_rows + ), f"Length of tracker angle list ({len(params.pv_array.tracker_angle)}) not equal to total number of PV tables ({params.pv_array.stream_rows * params.pv_array.span_rows})." + + tracker_angle_rad = np.radians(params.pv_array.tracker_angle[panel_id]) + else: + tracker_angle_rad = np.radians(params.pv_array.tracker_angle) + + if np.linalg.norm(phi) == 0: + K = np.zeros(params.pv_array.modules_per_span) + else: + K = np.abs( + 12 + * (np.delete(R_reaction_torque, params.pv_array.fixed_location)) + / ( + ( + params.pv_array.block_chord_div_by_panel_chord + * params.pv_array.panel_chord + ) + ** 3 + ) + / np.cos(tracker_angle_rad + phi) + / (np.sin(tracker_angle_rad + phi) - np.sin(tracker_angle_rad)) + / ( + params.pv_array.block_chord_div_by_panel_chord + * params.pv_array.panel_chord + / 2 + ) + ) + + # assume K is 0 at the fixed connector + K = np.flip(np.insert(K, params.pv_array.fixed_location, 0)) + + # K[0] is the stiffness at the most back connector (highest y) + # K[10] is the stiffness at the most front connector (lowest y) + # for block_bot_surface, it is numbered from lowest y to highest y, so flip K + + # if there are 10 modules per array, there are 11 connectors, K has shape of 10, the K at the fixed connector is not included. + for i in range(params.pv_array.modules_per_span + 1): + + name_K = f"spring_stiffness_{panel_id:.0f}_{i:.0f}" + setattr(self, name_K, K[i]) + + # num_panel_right_fixed = params.pv_array.modules_per_span // 2 + # num_panel_left_fixed = params.pv_array.modules_per_span - num_panel_right_fixed + + # # for right fixed part: + + # for panel_id in range(total_num_panels): + # total_torque_right_fixed = 0 + # total_torque_left_fixed = 0 + + # for i in range(num_panel_left_fixed, params.pv_array.modules_per_span): + # name = f"total_torque_panel_{panel_id:.0f}_{i:.0f}" + # total_torque_right_fixed += getattr(flow, name) + + # for i in range(num_panel_left_fixed): + # name = f"total_torque_panel_{panel_id:.0f}_{i:.0f}" + # total_torque_left_fixed += getattr(flow, name) + + # theo_matrix_right_fixed = np.zeros((num_panel_right_fixed+1, num_panel_right_fixed+1)) + # theo_vector_right_fixed = np.zeros(num_panel_right_fixed+1) + + # theo_matrix_right_fixed[0, :] = 1.0 + # theo_vector_right_fixed[0] = total_torque_right_fixed + # theo_matrix_right_fixed[1:, :-1] = params.pv_array.panel_span/params.pv_array.modules_per_span/Gt/Ipt + + # # as the double integral of torque along the span from front to back, left fixed part to right fixed part, it need to be flipped + # T_double_integral_right_fixed = [] + # T_right_fixed = [] + # for i in range(params.pv_array.modules_per_span): + # name = f"double_integral_total_torque_panel_{panel_id:.0f}_{params.pv_array.modules_per_span-1-i:.0f}" + # T_double_integral_right_fixed.append(getattr(flow, name)) + # name = f"total_torque_panel_{panel_id:.0f}_{params.pv_array.modules_per_span-1-i:.0f}" + # T_right_fixed.append(getattr(flow, name)) + # for i in range(num_panel_right_fixed): + # for j in range(i+1): + # theo_vector_right_fixed[i+1] += T_double_integral_right_fixed[-1-j]/(Gp*Ipp) + # theo_vector_right_fixed[i+1] -= (i+1-j)*T_right_fixed[-1-j]/(Gp*Ipp)*params.pv_array.panel_span/params.pv_array.modules_per_span + # for j in range(i): + # theo_matrix_right_fixed[i+1, :-1-j-1] += params.pv_array.panel_span/params.pv_array.modules_per_span/(Gt*Ipt) + + # for i in range(num_panel_right_fixed): + # for j in range(i+1): + # theo_matrix_right_fixed[i+1, num_panel_right_fixed-j:] -= params.pv_array.panel_span/params.pv_array.modules_per_span/Gp/Ipp + + # R_reaction_torque = np.dot(np.linalg.inv(theo_matrix_right_fixed), theo_vector_right_fixed) + + # tube_rotate_matrix = np.ones((num_panel_right_fixed, num_panel_right_fixed))*params.pv_array.panel_span/params.pv_array.modules_per_span*(1.0/Gt/Ipt) + + # for i in range(num_panel_right_fixed): + # for j in range(i): + # tube_rotate_matrix[i, :num_panel_right_fixed-1-j] += params.pv_array.panel_span/params.pv_array.modules_per_span*(1.0/Gt/Ipt) + + # # check the standalone code, why it need to be flipped. + # phi = np.flip(np.dot(tube_rotate_matrix, R_reaction_torque[:-1])) + + # # rotation of connector + # x_panel = np.arange(num_panel_right_fixed)*(params.pv_array.panel_span/params.pv_array.modules_per_span) # location of panels + + # block_length = params.pv_array.block_chord_div_by_panel_chord * params.pv_array.panel_chord + # block_width = params.pv_array.block_chord_div_by_panel_chord*params.pv_array.panel_span/params.pv_array.modules_per_span/2 + + # # To Do: check the angle, is this correct? + # array_rotation = (params.fluid.wind_direction + 90.0) % 360.0 + # array_rotation_rad = np.radians(array_rotation) + + # for i in range(num_panel_right_fixed): + + # phi_i = phi[i] + + # K_i = 12*(R_reaction_torque[i])/((block_length)**3)/np.cos(array_rotation_rad+phi_i)/(np.sin(array_rotation_rad+phi_i)-np.sin(array_rotation_rad))/block_width + + # name_K = f"spring_stiffness_{panel_id:.0f}_{params.pv_array.modules_per_span+1-i:.0f}" + # setattr(self, name_K, K_i) + + # # for left fixed part: + # theo_matrix_left_fixed = np.zeros((num_panel_left_fixed+1, num_panel_left_fixed+1)) + # theo_vector_left_fixed = np.zeros(num_panel_left_fixed+1) + # theo_matrix_left_fixed[0, :] = 1.0 + # theo_vector_left_fixed[0] = total_torque_left_fixed + + # T_double_integral_left_fixed = [] + # T_left_fixed = [] + + # for i in range(num_panel_left_fixed): + # name = f"double_integral_total_torque_panel_{panel_id:.0f}_{num_panel_left_fixed-1-i:.0f}" + # T_double_integral_left_fixed.append(getattr(flow, name)) + # name = f"total_torque_panel_{panel_id:.0f}_{num_panel_left_fixed-1-i:.0f}" + # T_left_fixed.append(getattr(flow, name)) + + # for i in range(num_panel_left_fixed): + # for j in range(i+1): + # theo_matrix_left_fixed[i+1, j] = (i+1-j)*params.pv_array.panel_span/params.pv_array.modules_per_span/(Gp*Ipp) + # theo_vector_left_fixed[i+1] += T_double_integral_left_fixed[j]/(Gp*Ipp) + # for j in range(i): + # theo_vector_left_fixed[i+1] += (i-j)*T_left_fixed[j]/(Gp*Ipp)*params.pv_array.panel_span/params.pv_array.modules_per_span + + # for i in range(num_panel_left_fixed): + # theo_matrix_left_fixed[i+1:, i+1:] -= params.pv_array.panel_span/params.pv_array.modules_per_span/Gt/Ipt + + # R_reaction_torque = np.dot(np.linalg.inv(theo_matrix_left_fixed), theo_vector_left_fixed) # this is the torque applied to tube + + # tube_rotate_matrix = np.zeros((num_panel_left_fixed, num_panel_left_fixed)) + # for i in range(num_panel_left_fixed): + # tube_rotate_matrix[i:, i:] += params.pv_array.panel_span/params.pv_array.modules_per_span/Gt/Ipt + + # phi = np.dot(tube_rotate_matrix, R_reaction_torque[1:]) + + # # rotation of connector + # x_panel = np.arange(num_panel_left_fixed)*(params.pv_array.panel_span/params.pv_array.modules_per_span)+params.pv_array.panel_span/params.pv_array.modules_per_span # location of panles + + # for i in range(num_panel_left_fixed): + + # phi_i = phi[i] + + # K_i = 12*(R_reaction_torque[i+1])/((block_length)**3)/np.cos(array_rotation_rad+phi_i)/(np.sin(array_rotation_rad+phi_i)-np.sin(array_rotation_rad))/block_width + + # name_K = f"spring_stiffness_{panel_id:.0f}_{num_panel_left_fixed-1-i:.0f}" + + # setattr(self, name_K, K_i) + def update_a(self, u, u_old, v_old, a_old, dt, beta, ufl=True): # Update formula for acceleration # a = 1/(2*beta)*((u - u0 - v0*dt)/(0.5*dt*dt) - (1-2*beta)*a0) @@ -126,7 +501,7 @@ def update_fields(self, u, u_old, v_old, a_old, dt, beta, gamma): def avg(self, x_old, x_new, alpha): return alpha * x_old + (1 - alpha) * x_new - def build_forms(self, domain, params, structure): + def build_forms(self, domain, params, structure, flow): """Builds all variational statements This method creates all the functions, expressions, and variational @@ -213,6 +588,9 @@ def c(u, u_): def k_nominal(u, u_): return ufl.inner(P_(u), ufl.grad(u_)) + def k_nominal_connector(u, u_): + return ufl.inner(P_connector(u), ufl.grad(u_)) + # The deformation gradient, F = I + dy/dX def F_(u): I = ufl.Identity(len(u)) @@ -242,6 +620,20 @@ def S_(u): S_svk = structure.lame_lambda * ufl.tr(E) * I + 2.0 * structure.lame_mu * E return S_svk + # The second Piola–Kirchhoff stress, S + def S_connector(u): + E = E_(u) + I = ufl.Identity(len(u)) + + # return lamda * ufl.tr(E) * I + 2.0 * mu * (E - ufl.tr(E) * I / 3.0) + # TODO: Why does the above form give a better result and where does it come from? + + S_svk = ( + structure.lame_lambda_connector * ufl.tr(E) * I + + 2.0 * structure.lame_mu_connector * E + ) + return S_svk + # The first Piola–Kirchhoff stress tensor, P = F*S def P_(u): F = F_(u) @@ -249,6 +641,12 @@ def P_(u): # return ufl.inv(F) * S return F * S + def P_connector(u): + F = F_(u) + S = S_connector(u) + # return ufl.inv(F) * S + return F * S + # self.uh_exp = dolfinx.fem.Function(self.V, name="Deformation") # def σ(v): @@ -302,14 +700,64 @@ def P_(u): F = ufl.grad(self.u) + ufl.Identity(len(self.u)) J = ufl.det(F) - self.res = ( - m(self.avg(self.a_old, a_new, self.alpha_m), self.u_) * ufl.dx - + c(self.avg(self.v_old, v_new, self.alpha_f), self.u_) * ufl.dx - + k_nominal(self.avg(self.u_old, self.u, self.alpha_f), self.u_) * ufl.dx - - structure.rho * ufl.inner(self.f, self.u_) * ufl.dx - - ufl.dot(ufl.dot(self.stress_predicted * J * ufl.inv(F.T), n), self.u_) - * self.ds - ) # - Wext(self.u) + + self.z_unit_vector = dolfinx.fem.Constant( + domain.structure.msh, [0.0, 0.0, 1.0] + ) # surface traction, N/m^2 + + if domain.modeling_torque_tube and params.general.geometry_module == "panels3d": + self.calculate_K_for_Robin_BC(domain, flow, params) + + dx_structure = ufl.Measure( + "dx", domain=domain.structure.msh, subdomain_data=domain.structure.cell_tags + ) + + if domain.modeling_torque_tube and params.general.geometry_module == "panels3d": + self.res = ( + m(self.avg(self.a_old, a_new, self.alpha_m), self.u_) * dx_structure + + c(self.avg(self.v_old, v_new, self.alpha_f), self.u_) * dx_structure + + k_nominal(self.avg(self.u_old, self.u, self.alpha_f), self.u_) + * dx_structure(domain.domain_markers["modules"]["idx"]) + + k_nominal_connector( + self.avg(self.u_old, self.u, self.alpha_f), self.u_ + ) + * dx_structure(domain.domain_markers["connectors"]["idx"]) + - structure.rho + * ufl.inner(self.f, self.u_) + * dx_structure(domain.domain_markers["modules"]["idx"]) + - structure.rho_connector + * ufl.inner(self.f, self.u_) + * dx_structure(domain.domain_markers["connectors"]["idx"]) + - ufl.dot(ufl.dot(self.stress_predicted * J * ufl.inv(F.T), n), self.u_) + * self.ds + ) # - Wext(self.u) + + # Robin boundary condition terms + for panel_id in range( + params.pv_array.stream_rows * params.pv_array.span_rows + ): + for i in range(params.pv_array.modules_per_span + 1): + name_K = f"spring_stiffness_{panel_id:.0f}_{i:.0f}" + K_springs = dolfinx.fem.Constant( + domain.structure.msh, float(getattr(self, name_K)) + ) + self.res -= ufl.dot( + K_springs * self.u_, self.z_unit_vector + ) * self.ds( + domain.domain_markers[f"block_bottom_{panel_id:.0f}_{i:.0f}"][ + "idx" + ] + ) + else: + self.res = ( + m(self.avg(self.a_old, a_new, self.alpha_m), self.u_) * dx_structure + + c(self.avg(self.v_old, v_new, self.alpha_f), self.u_) * dx_structure + + k_nominal(self.avg(self.u_old, self.u, self.alpha_f), self.u_) + * dx_structure + - structure.rho * ufl.inner(self.f, self.u_) * dx_structure + - ufl.dot(ufl.dot(self.stress_predicted * J * ufl.inv(F.T), n), self.u_) + * self.ds + ) # - Wext(self.u) # self.a = dolfinx.fem.form(ufl.lhs(res)) # self.L = dolfinx.fem.form(ufl.rhs(res)) @@ -484,26 +932,52 @@ def solve(self, params, dataIO, structure): try: idx = structure.north_east_corner_dofs[0] - # idx = self.north_east_corner_dofs[0] - nw_corner_accel = self.u.x.array[ + north_east_corner_deformation = self.u.x.array[ structure.ndim * idx : structure.ndim * idx + structure.ndim ].astype(np.float64) - print(nw_corner_accel) + print(f"Deformation: {north_east_corner_deformation}") + + except: + north_east_corner_deformation = np.zeros(structure.ndim, dtype=np.float64) + + try: + idx = structure.north_east_corner_dofs[0] + north_east_corner_acceleration = self.a_old.x.array[ + structure.ndim * idx : structure.ndim * idx + structure.ndim + ].astype(np.float64) + print(f"Acceleration: {north_east_corner_acceleration}") + except: - nw_corner_accel = np.zeros(structure.ndim, dtype=np.float64) + north_east_corner_acceleration = np.zeros(structure.ndim, dtype=np.float64) - nw_corner_accel_global = np.zeros( + # Initialize a buffer to collect everything into + north_east_corner_deformation_global = np.zeros( (self.num_procs, structure.ndim), dtype=np.float64 ) - self.comm.Gather(nw_corner_accel, nw_corner_accel_global, root=0) + north_east_corner_acceleration_global = np.zeros( + (self.num_procs, structure.ndim), dtype=np.float64 + ) - # print(f"Acceleration at North West corner = {nw_corner_accel}") + # Gather all points (many of which are zeros) to rank 0 + self.comm.Gather( + north_east_corner_deformation, north_east_corner_deformation_global, root=0 + ) + + self.comm.Gather( + north_east_corner_acceleration, + north_east_corner_acceleration_global, + root=0, + ) if self.rank == 0: - norm2 = np.sum(nw_corner_accel_global**2, axis=1) + norm2 = np.sum(north_east_corner_deformation_global**2, axis=1) + max_norm2_idx = np.argmax(norm2) + np_deformation = north_east_corner_deformation_global[max_norm2_idx, :] + + norm2 = np.sum(north_east_corner_acceleration_global**2, axis=1) max_norm2_idx = np.argmax(norm2) - np_accel = nw_corner_accel_global[max_norm2_idx, :] + np_acceleration = north_east_corner_acceleration_global[max_norm2_idx, :] accel_pos_filename = os.path.join( params.general.output_dir_sol, "accel_pos.csv" @@ -512,18 +986,35 @@ def solve(self, params, dataIO, structure): if self.first_call_to_solver: with open(accel_pos_filename, "w") as fp: - fp.write("#x-pos,y-pos,z-pos\n") if structure.ndim == 3: - fp.write(f"{np_accel[0]},{np_accel[1]},{np_accel[2]}\n") + fp.write( + "#x-deformation,y-deformation,z-deformation,x-acceleration,y-acceleration,z-acceleration\n" + ) + fp.write( + f"{np_deformation[0]},{np_deformation[1]},{np_deformation[2]}," + ) + fp.write( + f"{np_acceleration[0]},{np_acceleration[1]},{np_acceleration[2]}\n" + ) elif structure.ndim == 2: - fp.write(f"{np_accel[0]},{np_accel[1]}\n") + fp.write( + "#x-deformation,y-deformation,x-acceleration,y-acceleration\n" + ) + fp.write(f"{np_deformation[0]},{np_deformation[1]},") + fp.write(f"{np_acceleration[0]},{np_acceleration[1]}\n") else: with open(accel_pos_filename, "a") as fp: if structure.ndim == 3: - fp.write(f"{np_accel[0]},{np_accel[1]},{np_accel[2]}\n") + fp.write( + f"{np_deformation[0]},{np_deformation[1]},{np_deformation[2]}," + ) + fp.write( + f"{np_acceleration[0]},{np_acceleration[1]},{np_acceleration[2]}\n" + ) elif structure.ndim == 2: - fp.write(f"{np_accel[0]},{np_accel[1]}\n") + fp.write(f"{np_deformation[0]},{np_deformation[1]},") + fp.write(f"{np_acceleration[0]},{np_acceleration[1]}\n") if self.first_call_to_solver: self.first_call_to_solver = False diff --git a/pvade/structure/StructureMain.py b/pvade/structure/StructureMain.py index 6e55265d..3ca6d2fa 100644 --- a/pvade/structure/StructureMain.py +++ b/pvade/structure/StructureMain.py @@ -80,6 +80,10 @@ def __init__(self, domain, params): domain.structure.msh, params.structure.rho ) # Constant(0.) + self.rho_connector = dolfinx.fem.Constant( + domain.structure.msh, params.structure.rho_tube + ) + # Define structural properties self.E = params.structure.elasticity_modulus # 1.0e9 self.poissons_ratio = params.structure.poissons_ratio # 0.3 @@ -90,6 +94,21 @@ def __init__(self, domain, params): / ((1.0 + self.poissons_ratio) * (1.0 - 2.0 * self.poissons_ratio)) ) + # we are assuming the connectors has same mechanical properties as the tubes + self.E_connector = params.structure.elasticity_modulus_tube # 1.0e9 + self.poissons_ratio_connector = params.structure.poissons_ratio_tube # 0.3 + self.lame_mu_connector = self.E_connector / ( + 2.0 * (1.0 + self.poissons_ratio_connector) + ) + self.lame_lambda_connector = ( + self.E_connector + * self.poissons_ratio_connector + / ( + (1.0 + self.poissons_ratio_connector) + * (1.0 - 2.0 * self.poissons_ratio_connector) + ) + ) + if self.rank == 0: print( f"mu = {self.lame_mu} lambda = {self.lame_lambda} E = {self.E} nu = {self.poissons_ratio} density = {self.rho.value}" @@ -110,9 +129,21 @@ def _north_east_corner(x): else: tracker_angle_rad = np.radians(params.pv_array.tracker_angle) - x1 = 0.5 * params.pv_array.panel_chord * np.cos(tracker_angle_rad) - x2 = 0.5 * params.pv_array.panel_thickness * np.sin(tracker_angle_rad) - corner = [x1 - x2, 0.5 * params.pv_array.panel_span] + """ + y + + ^ + | + |---o NE corner (measured on bottom of panel) + | | + | | + |-----> x + """ + + corner = [ + 0.5 * params.pv_array.panel_chord * np.cos(tracker_angle_rad), + 0.5 * params.pv_array.panel_span, + ] east_edge = np.logical_and(corner[0] - eps < x[0], x[0] < corner[0] + eps) north_edge = np.logical_and(corner[1] - eps < x[1], x[1] < corner[1] + eps) @@ -121,12 +152,12 @@ def _north_east_corner(x): return north_east_corner - north_east_corner_facets = dolfinx.mesh.locate_entities_boundary( + north_east_corner_vertices = dolfinx.mesh.locate_entities_boundary( domain.structure.msh, 0, _north_east_corner ) self.north_east_corner_dofs = dolfinx.fem.locate_dofs_topological( - self.elasticity.V, 0, north_east_corner_facets + self.elasticity.V, 0, north_east_corner_vertices ) def build_boundary_conditions(self, domain, params): @@ -189,7 +220,7 @@ def update_fields(self, u, u_old, v_old, a_old, dt, beta, gamma): def avg(self, x_old, x_new, alpha): return alpha * x_old + (1 - alpha) * x_new - def build_forms(self, domain, params): + def build_forms(self, domain, params, flow): """Builds all variational statements This method creates all the functions, expressions, and variational @@ -205,7 +236,7 @@ def build_forms(self, domain, params): params (:obj:`pvade.Parameters.SimParams`): A SimParams object """ - self.elasticity.build_forms(domain, params, self) + self.elasticity.build_forms(domain, params, self, flow) def _assemble_system(self, params): """Pre-assemble all LHS matrices and RHS vectors diff --git a/pvade/structure/boundary_conditions.py b/pvade/structure/boundary_conditions.py index 0aa0ea0e..57fc98d1 100644 --- a/pvade/structure/boundary_conditions.py +++ b/pvade/structure/boundary_conditions.py @@ -272,8 +272,11 @@ def build_structure_boundary_conditions(domain, params, functionspace): total_num_panels = params.pv_array.stream_rows * params.pv_array.span_rows for num_panel in range(total_num_panels): - for location in params.structure.bc_list: - location_panel = f"{location}_{num_panel}" + for location in params.structure.bc_list: # it is empty + if location == "panel_top" or location == "panel_bottom": + location_panel = f"{location}_{num_panel}_0" + else: + location_panel = f"{location}_{num_panel}" # f"front_{num_panel}" , f"back_{num_panel}": # for location in [f"left_{num_panel}"]:# , f"right_{num_panel}": # for location in f"left_{num_panel}": @@ -344,43 +347,78 @@ def connection_point_up_helper(nodes_to_pin_between): return fn_handle - # Start pinning along the lines expressed byt numpy_pt_total_array - # First, determine the total number of rows (number of lines to pin) - num_nodes = np.shape(domain.numpy_pt_total_array)[0] + if domain.modeling_torque_tube and params.general.geometry_module == "panels3d": + # # Start pinning along the lines expressed byt numpy_pt_total_array + # The center line of connectors bottom surface is fixed to remove rigid body motion - # Determine how many pinning lines exist per each panel (e.g., 24 lines distributed on 8 panels means 3 lines per panel) - nodes_per_panel = int(num_nodes / total_num_panels) + if params.structure.tube_connection == True: + # If making torque tube connections, pass only those pinning lines to the BC identification function + tube_nodes = domain.numpy_pt_total_array[:, :] - # The torque tube entry (oriented spanwise along the middle, divides panel into upstream and downstream rectangular halves) - # is always the first entry, e.g., [0, ..., ..., 3, ..., ..., 6, ..., ...], [0, 3, 6] are the torque tubes - tube_nodes_idx = np.arange(0, num_nodes, nodes_per_panel, dtype=np.int64) + facet_uppoint = dolfinx.mesh.locate_entities( + domain.structure.msh, 1, connection_point_up_helper(tube_nodes) + ) + dofs_disp = dolfinx.fem.locate_dofs_topological( + functionspace, 1, [facet_uppoint] + ) - if params.structure.tube_connection == True: - # If making torque tube connections, pass only those pinning lines to the BC identification function - tube_nodes = domain.numpy_pt_total_array[tube_nodes_idx, :] + bc.append(dolfinx.fem.dirichletbc(zero_vec, dofs_disp, functionspace)) - facet_uppoint = dolfinx.mesh.locate_entities( - domain.structure.msh, 1, connection_point_up_helper(tube_nodes) - ) - dofs_disp = dolfinx.fem.locate_dofs_topological( - functionspace, 1, [facet_uppoint] - ) + if params.structure.motor_connection == True: + # The bottom surface of the center connector is fixed to represent the motor mount + motor_location = ( + params.pv_array.fixed_location + ) # fixed location along the span, if fixed_location=5, it means the left boundary of the 6th connector is fixed + for panel_id in range(total_num_panels): - bc.append(dolfinx.fem.dirichletbc(zero_vec, dofs_disp, functionspace)) + mount_facet = domain.structure.facet_tags.find( + domain.domain_markers[ + f"block_left_{panel_id:.0f}_{motor_location:.0f}" + ]["idx"] + ) - if params.structure.motor_connection == True: - # If making motor mount connections, pass only those pinning lines to the BC identification function - # this is done by making a copy of the numpy_pt_total_array with the torque tube lines *deleted* - # not done in place, so numpy_pt_total_array remains unaltered. - motor_nodes = np.delete(domain.numpy_pt_total_array, tube_nodes_idx, axis=0) + dofs_disp = dolfinx.fem.locate_dofs_topological( + functionspace, 2, [mount_facet] + ) - facet_uppoint = dolfinx.mesh.locate_entities( - domain.structure.msh, 1, connection_point_up_helper(motor_nodes) - ) - dofs_disp = dolfinx.fem.locate_dofs_topological( - functionspace, 1, [facet_uppoint] - ) + bc.append(dolfinx.fem.dirichletbc(zero_vec, dofs_disp, functionspace)) + else: + # Start pinning along the lines expressed byt numpy_pt_total_array + # First, determine the total number of rows (number of lines to pin) + num_nodes = np.shape(domain.numpy_pt_total_array)[0] + + # Determine how many pinning lines exist per each panel (e.g., 24 lines distributed on 8 panels means 3 lines per panel) + nodes_per_panel = int(num_nodes / total_num_panels) - bc.append(dolfinx.fem.dirichletbc(zero_vec, dofs_disp, functionspace)) + # The torque tube entry (oriented spanwise along the middle, divides panel into upstream and downstream rectangular halves) + # is always the first entry, e.g., [0, ..., ..., 3, ..., ..., 6, ..., ...], [0, 3, 6] are the torque tubes + tube_nodes_idx = np.arange(0, num_nodes, nodes_per_panel, dtype=np.int64) + + if params.structure.tube_connection == True: + # If making torque tube connections, pass only those pinning lines to the BC identification function + tube_nodes = domain.numpy_pt_total_array[tube_nodes_idx, :] + + facet_uppoint = dolfinx.mesh.locate_entities( + domain.structure.msh, 1, connection_point_up_helper(tube_nodes) + ) + dofs_disp = dolfinx.fem.locate_dofs_topological( + functionspace, 1, [facet_uppoint] + ) + + bc.append(dolfinx.fem.dirichletbc(zero_vec, dofs_disp, functionspace)) + + if params.structure.motor_connection == True: + # If making motor mount connections, pass only those pinning lines to the BC identification function + # this is done by making a copy of the numpy_pt_total_array with the torque tube lines *deleted* + # not done in place, so numpy_pt_total_array remains unaltered. + motor_nodes = np.delete(domain.numpy_pt_total_array, tube_nodes_idx, axis=0) + + facet_uppoint = dolfinx.mesh.locate_entities( + domain.structure.msh, 1, connection_point_up_helper(motor_nodes) + ) + dofs_disp = dolfinx.fem.locate_dofs_topological( + functionspace, 1, [facet_uppoint] + ) + bc.append(dolfinx.fem.dirichletbc(zero_vec, dofs_disp, functionspace)) return bc diff --git a/pvade/tests/input/mesh/panels3d/domain_markers.yaml b/pvade/tests/input/mesh/panels3d/domain_markers.yaml index 2af5e530..50b75e6d 100644 --- a/pvade/tests/input/mesh/panels3d/domain_markers.yaml +++ b/pvade/tests/input/mesh/panels3d/domain_markers.yaml @@ -1,110 +1,56 @@ -_current_idx: 135 -back_0: - entity: facet +_current_idx: 51 +fluid: + entity: cell gmsh_tags: - 8 - idx: 2 -back_1: - entity: facet - gmsh_tags: - - 17 - idx: 8 -back_10: - entity: facet - gmsh_tags: - - 98 - idx: 62 -back_11: - entity: facet - gmsh_tags: - - 107 - idx: 68 -back_12: - entity: facet - gmsh_tags: - - 116 - idx: 74 -back_13: - entity: facet - gmsh_tags: - - 125 - idx: 80 -back_14: - entity: facet - gmsh_tags: - - 134 - idx: 86 -back_15: - entity: facet - gmsh_tags: - - 143 - idx: 92 -back_16: - entity: facet - gmsh_tags: - - 152 - idx: 98 -back_17: - entity: facet - gmsh_tags: - - 161 - idx: 104 -back_18: - entity: facet - gmsh_tags: - - 170 - idx: 110 -back_19: - entity: facet - gmsh_tags: - - 179 - idx: 116 -back_2: - entity: facet - gmsh_tags: - - 26 - idx: 14 -back_20: - entity: facet + idx: 50 +modules: + entity: cell gmsh_tags: - - 188 - idx: 122 -back_3: + - 1 + - 2 + - 3 + - 4 + - 5 + - 6 + - 7 + idx: 49 +panel_back_0: entity: facet gmsh_tags: - - 35 - idx: 20 -back_4: + - 10 + idx: 4 +panel_back_1: entity: facet gmsh_tags: - - 44 - idx: 26 -back_5: + - 19 + idx: 10 +panel_back_2: entity: facet gmsh_tags: - - 53 - idx: 32 -back_6: + - 28 + idx: 16 +panel_back_3: entity: facet gmsh_tags: - - 62 - idx: 38 -back_7: + - 37 + idx: 22 +panel_back_4: entity: facet gmsh_tags: - - 71 - idx: 44 -back_8: + - 46 + idx: 28 +panel_back_5: entity: facet gmsh_tags: - - 80 - idx: 50 -back_9: + - 55 + idx: 34 +panel_back_6: entity: facet gmsh_tags: - - 89 - idx: 56 -bottom_0: + - 64 + idx: 40 +panel_bottom_0_0: entity: facet gmsh_tags: - 11 @@ -112,7 +58,7 @@ bottom_0: - 13 - 14 idx: 5 -bottom_1: +panel_bottom_1_0: entity: facet gmsh_tags: - 20 @@ -120,87 +66,7 @@ bottom_1: - 22 - 23 idx: 11 -bottom_10: - entity: facet - gmsh_tags: - - 101 - - 102 - - 103 - - 104 - idx: 65 -bottom_11: - entity: facet - gmsh_tags: - - 110 - - 111 - - 112 - - 113 - idx: 71 -bottom_12: - entity: facet - gmsh_tags: - - 119 - - 120 - - 121 - - 122 - idx: 77 -bottom_13: - entity: facet - gmsh_tags: - - 128 - - 129 - - 130 - - 131 - idx: 83 -bottom_14: - entity: facet - gmsh_tags: - - 137 - - 138 - - 139 - - 140 - idx: 89 -bottom_15: - entity: facet - gmsh_tags: - - 146 - - 147 - - 148 - - 149 - idx: 95 -bottom_16: - entity: facet - gmsh_tags: - - 155 - - 156 - - 157 - - 158 - idx: 101 -bottom_17: - entity: facet - gmsh_tags: - - 164 - - 165 - - 166 - - 167 - idx: 107 -bottom_18: - entity: facet - gmsh_tags: - - 173 - - 174 - - 175 - - 176 - idx: 113 -bottom_19: - entity: facet - gmsh_tags: - - 182 - - 183 - - 184 - - 185 - idx: 119 -bottom_2: +panel_bottom_2_0: entity: facet gmsh_tags: - 29 @@ -208,15 +74,7 @@ bottom_2: - 31 - 32 idx: 17 -bottom_20: - entity: facet - gmsh_tags: - - 191 - - 192 - - 193 - - 194 - idx: 125 -bottom_3: +panel_bottom_3_0: entity: facet gmsh_tags: - 38 @@ -224,7 +82,7 @@ bottom_3: - 40 - 41 idx: 23 -bottom_4: +panel_bottom_4_0: entity: facet gmsh_tags: - 47 @@ -232,7 +90,7 @@ bottom_4: - 49 - 50 idx: 29 -bottom_5: +panel_bottom_5_0: entity: facet gmsh_tags: - 56 @@ -240,7 +98,7 @@ bottom_5: - 58 - 59 idx: 35 -bottom_6: +panel_bottom_6_0: entity: facet gmsh_tags: - 65 @@ -248,507 +106,173 @@ bottom_6: - 67 - 68 idx: 41 -bottom_7: - entity: facet - gmsh_tags: - - 74 - - 75 - - 76 - - 77 - idx: 47 -bottom_8: - entity: facet - gmsh_tags: - - 83 - - 84 - - 85 - - 86 - idx: 53 -bottom_9: - entity: facet - gmsh_tags: - - 92 - - 93 - - 94 - - 95 - idx: 59 -fluid: - entity: cell - gmsh_tags: - - 22 - idx: 134 -front_0: - entity: facet - gmsh_tags: - - 7 - idx: 1 -front_1: - entity: facet - gmsh_tags: - - 16 - idx: 7 -front_10: - entity: facet - gmsh_tags: - - 97 - idx: 61 -front_11: - entity: facet - gmsh_tags: - - 106 - idx: 67 -front_12: - entity: facet - gmsh_tags: - - 115 - idx: 73 -front_13: - entity: facet - gmsh_tags: - - 124 - idx: 79 -front_14: - entity: facet - gmsh_tags: - - 133 - idx: 85 -front_15: - entity: facet - gmsh_tags: - - 142 - idx: 91 -front_16: - entity: facet - gmsh_tags: - - 151 - idx: 97 -front_17: - entity: facet - gmsh_tags: - - 160 - idx: 103 -front_18: - entity: facet - gmsh_tags: - - 169 - idx: 109 -front_19: - entity: facet - gmsh_tags: - - 178 - idx: 115 -front_2: - entity: facet - gmsh_tags: - - 25 - idx: 13 -front_20: - entity: facet - gmsh_tags: - - 187 - idx: 121 -front_3: - entity: facet - gmsh_tags: - - 34 - idx: 19 -front_4: - entity: facet - gmsh_tags: - - 43 - idx: 25 -front_5: - entity: facet - gmsh_tags: - - 52 - idx: 31 -front_6: - entity: facet - gmsh_tags: - - 61 - idx: 37 -front_7: - entity: facet - gmsh_tags: - - 70 - idx: 43 -front_8: - entity: facet - gmsh_tags: - - 79 - idx: 49 -front_9: - entity: facet - gmsh_tags: - - 88 - idx: 55 -left_0: +panel_front_0: entity: facet gmsh_tags: - 9 idx: 3 -left_1: +panel_front_1: entity: facet gmsh_tags: - 18 idx: 9 -left_10: - entity: facet - gmsh_tags: - - 99 - idx: 63 -left_11: - entity: facet - gmsh_tags: - - 108 - idx: 69 -left_12: - entity: facet - gmsh_tags: - - 117 - idx: 75 -left_13: - entity: facet - gmsh_tags: - - 126 - idx: 81 -left_14: - entity: facet - gmsh_tags: - - 135 - idx: 87 -left_15: - entity: facet - gmsh_tags: - - 144 - idx: 93 -left_16: - entity: facet - gmsh_tags: - - 153 - idx: 99 -left_17: - entity: facet - gmsh_tags: - - 162 - idx: 105 -left_18: - entity: facet - gmsh_tags: - - 171 - idx: 111 -left_19: - entity: facet - gmsh_tags: - - 180 - idx: 117 -left_2: +panel_front_2: entity: facet gmsh_tags: - 27 idx: 15 -left_20: - entity: facet - gmsh_tags: - - 189 - idx: 123 -left_3: +panel_front_3: entity: facet gmsh_tags: - 36 idx: 21 -left_4: +panel_front_4: entity: facet gmsh_tags: - 45 idx: 27 -left_5: +panel_front_5: entity: facet gmsh_tags: - 54 idx: 33 -left_6: +panel_front_6: entity: facet gmsh_tags: - 63 idx: 39 -left_7: - entity: facet - gmsh_tags: - - 72 - idx: 45 -left_8: - entity: facet - gmsh_tags: - - 81 - idx: 51 -left_9: - entity: facet - gmsh_tags: - - 90 - idx: 57 -right_0: - entity: facet - gmsh_tags: - - 10 - idx: 4 -right_1: - entity: facet - gmsh_tags: - - 19 - idx: 10 -right_10: +panel_left_0: entity: facet gmsh_tags: - - 100 - idx: 64 -right_11: - entity: facet - gmsh_tags: - - 109 - idx: 70 -right_12: - entity: facet - gmsh_tags: - - 118 - idx: 76 -right_13: - entity: facet - gmsh_tags: - - 127 - idx: 82 -right_14: - entity: facet - gmsh_tags: - - 136 - idx: 88 -right_15: - entity: facet - gmsh_tags: - - 145 - idx: 94 -right_16: + - 7 + idx: 1 +panel_left_1: entity: facet gmsh_tags: - - 154 - idx: 100 -right_17: + - 16 + idx: 7 +panel_left_2: entity: facet gmsh_tags: - - 163 - idx: 106 -right_18: + - 25 + idx: 13 +panel_left_3: entity: facet gmsh_tags: - - 172 - idx: 112 -right_19: + - 34 + idx: 19 +panel_left_4: entity: facet gmsh_tags: - - 181 - idx: 118 -right_2: + - 43 + idx: 25 +panel_left_5: entity: facet gmsh_tags: - - 28 - idx: 16 -right_20: + - 52 + idx: 31 +panel_left_6: entity: facet gmsh_tags: - - 190 - idx: 124 -right_3: + - 61 + idx: 37 +panel_right_0: entity: facet gmsh_tags: - - 37 - idx: 22 -right_4: + - 8 + idx: 2 +panel_right_1: entity: facet gmsh_tags: - - 46 - idx: 28 -right_5: + - 17 + idx: 8 +panel_right_2: entity: facet gmsh_tags: - - 55 - idx: 34 -right_6: + - 26 + idx: 14 +panel_right_3: entity: facet gmsh_tags: - - 64 - idx: 40 -right_7: + - 35 + idx: 20 +panel_right_4: entity: facet gmsh_tags: - - 73 - idx: 46 -right_8: + - 44 + idx: 26 +panel_right_5: entity: facet gmsh_tags: - - 82 - idx: 52 -right_9: + - 53 + idx: 32 +panel_right_6: entity: facet gmsh_tags: - - 91 - idx: 58 -structure: - entity: cell - gmsh_tags: - - 1 - - 2 - - 3 - - 4 - - 5 - - 6 - - 7 - - 8 - - 9 - - 10 - - 11 - - 12 - - 13 - - 14 - - 15 - - 16 - - 17 - - 18 - - 19 - - 20 - - 21 - idx: 133 -top_0: + - 62 + idx: 38 +panel_top_0_0: entity: facet gmsh_tags: - 15 idx: 6 -top_1: +panel_top_1_0: entity: facet gmsh_tags: - 24 idx: 12 -top_10: - entity: facet - gmsh_tags: - - 105 - idx: 66 -top_11: - entity: facet - gmsh_tags: - - 114 - idx: 72 -top_12: - entity: facet - gmsh_tags: - - 123 - idx: 78 -top_13: - entity: facet - gmsh_tags: - - 132 - idx: 84 -top_14: - entity: facet - gmsh_tags: - - 141 - idx: 90 -top_15: - entity: facet - gmsh_tags: - - 150 - idx: 96 -top_16: - entity: facet - gmsh_tags: - - 159 - idx: 102 -top_17: - entity: facet - gmsh_tags: - - 168 - idx: 108 -top_18: - entity: facet - gmsh_tags: - - 177 - idx: 114 -top_19: - entity: facet - gmsh_tags: - - 186 - idx: 120 -top_2: +panel_top_2_0: entity: facet gmsh_tags: - 33 idx: 18 -top_20: - entity: facet - gmsh_tags: - - 195 - idx: 126 -top_3: +panel_top_3_0: entity: facet gmsh_tags: - 42 idx: 24 -top_4: +panel_top_4_0: entity: facet gmsh_tags: - 51 idx: 30 -top_5: +panel_top_5_0: entity: facet gmsh_tags: - 60 idx: 36 -top_6: +panel_top_6_0: entity: facet gmsh_tags: - 69 idx: 42 -top_7: - entity: facet - gmsh_tags: - - 78 - idx: 48 -top_8: - entity: facet - gmsh_tags: - - 87 - idx: 54 -top_9: - entity: facet - gmsh_tags: - - 96 - idx: 60 x_max: entity: facet gmsh_tags: - - 6 - idx: 132 + - 75 + idx: 48 x_min: entity: facet gmsh_tags: - - 1 - idx: 127 + - 70 + idx: 43 y_max: entity: facet gmsh_tags: - - 4 - idx: 130 + - 73 + idx: 46 y_min: entity: facet gmsh_tags: - - 2 - idx: 128 + - 71 + idx: 44 z_max: entity: facet gmsh_tags: - - 3 - idx: 129 + - 72 + idx: 45 z_min: entity: facet gmsh_tags: - - 5 - idx: 131 + - 74 + idx: 47 diff --git a/pvade/tests/input/mesh/panels3d/fluid_mesh.h5 b/pvade/tests/input/mesh/panels3d/fluid_mesh.h5 index d382a0c5..4f3f1571 100644 Binary files a/pvade/tests/input/mesh/panels3d/fluid_mesh.h5 and b/pvade/tests/input/mesh/panels3d/fluid_mesh.h5 differ diff --git a/pvade/tests/input/mesh/panels3d/fluid_mesh.xdmf b/pvade/tests/input/mesh/panels3d/fluid_mesh.xdmf index c60059dc..ba06796e 100644 --- a/pvade/tests/input/mesh/panels3d/fluid_mesh.xdmf +++ b/pvade/tests/input/mesh/panels3d/fluid_mesh.xdmf @@ -3,29 +3,29 @@ - - fluid_mesh.h5:/Mesh/fluid_mesh.xdmf/topology + + fluid_mesh.h5:/Mesh/fluid_mesh.xdmf/topology - fluid_mesh.h5:/Mesh/fluid_mesh.xdmf/geometry + fluid_mesh.h5:/Mesh/fluid_mesh.xdmf/geometry - - fluid_mesh.h5:/MeshTags/cell_tags/topology + + fluid_mesh.h5:/MeshTags/cell_tags/topology - fluid_mesh.h5:/MeshTags/cell_tags/Values + fluid_mesh.h5:/MeshTags/cell_tags/Values - - fluid_mesh.h5:/MeshTags/facet_tags/topology + + fluid_mesh.h5:/MeshTags/facet_tags/topology - fluid_mesh.h5:/MeshTags/facet_tags/Values + fluid_mesh.h5:/MeshTags/facet_tags/Values diff --git a/pvade/tests/input/mesh/panels3d/numpy_fixation_points.csv b/pvade/tests/input/mesh/panels3d/numpy_fixation_points.csv index 50ae95ab..cc278eee 100644 --- a/pvade/tests/input/mesh/panels3d/numpy_fixation_points.csv +++ b/pvade/tests/input/mesh/panels3d/numpy_fixation_points.csv @@ -1,18 +1,4 @@ # start_x,start_y,start_z,end_x,end_y,end_z, -2.499999999999857891e-02,-1.050000000000000000e+01,1.456698729810778081e+00,2.499999999999857891e-02,-3.500000000000000000e+00,1.456698729810778081e+00 --8.410254037844389075e-01,-7.000000000000000000e+00,9.566987298107780813e-01,8.910254037844396180e-01,-7.000000000000000000e+00,1.956698729810778081e+00 -7.025000000000000355e+00,-1.050000000000000000e+01,1.456698729810778081e+00,7.025000000000000355e+00,-3.500000000000000000e+00,1.456698729810778081e+00 -6.158974596215561093e+00,-7.000000000000000000e+00,9.566987298107780813e-01,7.891025403784439618e+00,-7.000000000000000000e+00,1.956698729810778081e+00 -1.402500000000000036e+01,-1.050000000000000000e+01,1.456698729810778081e+00,1.402500000000000036e+01,-3.500000000000000000e+00,1.456698729810778081e+00 -1.315897459621556109e+01,-7.000000000000000000e+00,9.566987298107780813e-01,1.489102540378443962e+01,-7.000000000000000000e+00,1.956698729810778081e+00 -2.102499999999999858e+01,-1.050000000000000000e+01,1.456698729810778081e+00,2.102499999999999858e+01,-3.500000000000000000e+00,1.456698729810778081e+00 -2.015897459621556109e+01,-7.000000000000000000e+00,9.566987298107780813e-01,2.189102540378443962e+01,-7.000000000000000000e+00,1.956698729810778081e+00 -2.802499999999999858e+01,-1.050000000000000000e+01,1.456698729810778081e+00,2.802499999999999858e+01,-3.500000000000000000e+00,1.456698729810778081e+00 -2.715897459621556109e+01,-7.000000000000000000e+00,9.566987298107780813e-01,2.889102540378443962e+01,-7.000000000000000000e+00,1.956698729810778081e+00 -3.502499999999999858e+01,-1.050000000000000000e+01,1.456698729810778081e+00,3.502499999999999858e+01,-3.500000000000000000e+00,1.456698729810778081e+00 -3.415897459621556465e+01,-7.000000000000000000e+00,9.566987298107780813e-01,3.589102540378443962e+01,-7.000000000000000000e+00,1.956698729810778081e+00 -4.202499999999999858e+01,-1.050000000000000000e+01,1.456698729810778081e+00,4.202499999999999858e+01,-3.500000000000000000e+00,1.456698729810778081e+00 -4.115897459621556465e+01,-7.000000000000000000e+00,9.566987298107780813e-01,4.289102540378443962e+01,-7.000000000000000000e+00,1.956698729810778081e+00 2.499999999999857891e-02,-3.500000000000000000e+00,1.456698729810778081e+00,2.499999999999857891e-02,3.500000000000000000e+00,1.456698729810778081e+00 -8.410254037844389075e-01,0.000000000000000000e+00,9.566987298107780813e-01,8.910254037844396180e-01,0.000000000000000000e+00,1.956698729810778081e+00 7.025000000000000355e+00,-3.500000000000000000e+00,1.456698729810778081e+00,7.025000000000000355e+00,3.500000000000000000e+00,1.456698729810778081e+00 @@ -27,17 +13,3 @@ 3.415897459621556465e+01,0.000000000000000000e+00,9.566987298107780813e-01,3.589102540378443962e+01,0.000000000000000000e+00,1.956698729810778081e+00 4.202499999999999858e+01,-3.500000000000000000e+00,1.456698729810778081e+00,4.202499999999999858e+01,3.500000000000000000e+00,1.456698729810778081e+00 4.115897459621556465e+01,0.000000000000000000e+00,9.566987298107780813e-01,4.289102540378443962e+01,0.000000000000000000e+00,1.956698729810778081e+00 -2.499999999999857891e-02,3.500000000000000000e+00,1.456698729810778081e+00,2.499999999999857891e-02,1.050000000000000000e+01,1.456698729810778081e+00 --8.410254037844389075e-01,7.000000000000000000e+00,9.566987298107780813e-01,8.910254037844396180e-01,7.000000000000000000e+00,1.956698729810778081e+00 -7.025000000000000355e+00,3.500000000000000000e+00,1.456698729810778081e+00,7.025000000000000355e+00,1.050000000000000000e+01,1.456698729810778081e+00 -6.158974596215561093e+00,7.000000000000000000e+00,9.566987298107780813e-01,7.891025403784439618e+00,7.000000000000000000e+00,1.956698729810778081e+00 -1.402500000000000036e+01,3.500000000000000000e+00,1.456698729810778081e+00,1.402500000000000036e+01,1.050000000000000000e+01,1.456698729810778081e+00 -1.315897459621556109e+01,7.000000000000000000e+00,9.566987298107780813e-01,1.489102540378443962e+01,7.000000000000000000e+00,1.956698729810778081e+00 -2.102499999999999858e+01,3.500000000000000000e+00,1.456698729810778081e+00,2.102499999999999858e+01,1.050000000000000000e+01,1.456698729810778081e+00 -2.015897459621556109e+01,7.000000000000000000e+00,9.566987298107780813e-01,2.189102540378443962e+01,7.000000000000000000e+00,1.956698729810778081e+00 -2.802499999999999858e+01,3.500000000000000000e+00,1.456698729810778081e+00,2.802499999999999858e+01,1.050000000000000000e+01,1.456698729810778081e+00 -2.715897459621556109e+01,7.000000000000000000e+00,9.566987298107780813e-01,2.889102540378443962e+01,7.000000000000000000e+00,1.956698729810778081e+00 -3.502499999999999858e+01,3.500000000000000000e+00,1.456698729810778081e+00,3.502499999999999858e+01,1.050000000000000000e+01,1.456698729810778081e+00 -3.415897459621556465e+01,7.000000000000000000e+00,9.566987298107780813e-01,3.589102540378443962e+01,7.000000000000000000e+00,1.956698729810778081e+00 -4.202499999999999858e+01,3.500000000000000000e+00,1.456698729810778081e+00,4.202499999999999858e+01,1.050000000000000000e+01,1.456698729810778081e+00 -4.115897459621556465e+01,7.000000000000000000e+00,9.566987298107780813e-01,4.289102540378443962e+01,7.000000000000000000e+00,1.956698729810778081e+00 diff --git a/pvade/tests/input/mesh/panels3d/structure_mesh.h5 b/pvade/tests/input/mesh/panels3d/structure_mesh.h5 index d4acd84e..9c47e339 100644 Binary files a/pvade/tests/input/mesh/panels3d/structure_mesh.h5 and b/pvade/tests/input/mesh/panels3d/structure_mesh.h5 differ diff --git a/pvade/tests/input/mesh/panels3d/structure_mesh.xdmf b/pvade/tests/input/mesh/panels3d/structure_mesh.xdmf index df8c6c3d..a0748306 100644 --- a/pvade/tests/input/mesh/panels3d/structure_mesh.xdmf +++ b/pvade/tests/input/mesh/panels3d/structure_mesh.xdmf @@ -3,29 +3,29 @@ - - structure_mesh.h5:/Mesh/structure_mesh.xdmf/topology + + structure_mesh.h5:/Mesh/structure_mesh.xdmf/topology - structure_mesh.h5:/Mesh/structure_mesh.xdmf/geometry + structure_mesh.h5:/Mesh/structure_mesh.xdmf/geometry - - structure_mesh.h5:/MeshTags/cell_tags/topology + + structure_mesh.h5:/MeshTags/cell_tags/topology - structure_mesh.h5:/MeshTags/cell_tags/Values + structure_mesh.h5:/MeshTags/cell_tags/Values - - structure_mesh.h5:/MeshTags/facet_tags/topology + + structure_mesh.h5:/MeshTags/facet_tags/topology - structure_mesh.h5:/MeshTags/facet_tags/Values + structure_mesh.h5:/MeshTags/facet_tags/Values diff --git a/pvade/tests/input/yaml/embedded_box.yaml b/pvade/tests/input/yaml/embedded_box.yaml index 0bd12d54..966a1473 100644 --- a/pvade/tests/input/yaml/embedded_box.yaml +++ b/pvade/tests/input/yaml/embedded_box.yaml @@ -9,18 +9,19 @@ domain: x_max: 1.0 y_min: -1.0 y_max: 1.0 - z_min: -1.0 - z_max: 1.0 + z_min: 0 + z_max: 2.0 l_char: 0.05 pv_array: stream_rows: 1 span_rows: 1 stream_spacing: 1.0 span_spacing: 1.0 - elevation: 0.0 + elevation: 1. panel_chord: 1.2 panel_span: 1.1 panel_thickness: 1.0 tracker_angle: 0.0 + # span_fixation_pts: null diff --git a/pvade/tests/input/yaml/flag2d.yaml b/pvade/tests/input/yaml/flag2d.yaml index 02124225..e1a8eae7 100644 --- a/pvade/tests/input/yaml/flag2d.yaml +++ b/pvade/tests/input/yaml/flag2d.yaml @@ -4,7 +4,7 @@ general: mesh_only: false structural_analysis: True fluid_analysis: True - input_mesh_dir: pvade/tests/input/mesh/flag2d + input_mesh_dir: #pvade/tests/input/mesh/flag2d domain: x_min: 0.0 x_max: 2.5 @@ -53,7 +53,7 @@ structure: body_force_x: 0 body_force_y: 0 body_force_z: 0 #100 - bc_list: ["left"] + bc_list: ["panel_left"] motor_connection: False tube_connection: False beta_relaxation: 0.005 diff --git a/pvade/tests/test_fsi_mesh.py b/pvade/tests/test_fsi_mesh.py index 2d4f4758..24c8ebd6 100644 --- a/pvade/tests/test_fsi_mesh.py +++ b/pvade/tests/test_fsi_mesh.py @@ -132,6 +132,8 @@ def test_meshing_3dpanels_rotations(wind_direction, num_stream_rows, num_span_ro # Arrays always start at xc = 0, but are centered in the y-direction # so now shift the mean to 0.0 + # yc -= np.mean(yc) + yc = yc.astype(float) yc -= np.mean(yc) counter = 0 diff --git a/pvade/tests/test_input_files.py b/pvade/tests/test_input_files.py index 01c13c86..77af7459 100644 --- a/pvade/tests/test_input_files.py +++ b/pvade/tests/test_input_files.py @@ -19,7 +19,7 @@ def launch_sim(test_file): tf = dt * 10 # ten timesteps command = ( - f"mpirun -n 8 python " + f"mpirun -n 4 python -u " + rootdir + "/pvade_main.py --input_file " + test_file["path_to_file"] diff --git a/pvade/tests/test_mesh_movement.py b/pvade/tests/test_mesh_movement.py index 4de08c88..f52daf5f 100644 --- a/pvade/tests/test_mesh_movement.py +++ b/pvade/tests/test_mesh_movement.py @@ -49,7 +49,11 @@ def test_calc_distance_to_panel_surface(): dx = params.domain.x_max - 0.5 * params.pv_array.panel_chord dy = params.domain.y_max - 0.5 * params.pv_array.panel_span - dz = params.domain.z_max - 0.5 * params.pv_array.panel_thickness + dz = ( + params.domain.z_max + - 0.5 * params.pv_array.panel_thickness + - params.pv_array.elevation + ) truth_max_dist = np.sqrt(dx * dx + dy * dy + dz * dz) assert np.isclose(max_dist, truth_max_dist) @@ -114,7 +118,8 @@ def __init__(self, domain, x_shift, y_shift, z_shift): np.amin(structure_coords_before[:, 1]), -0.5 * params.pv_array.panel_span ) assert np.isclose( - np.amin(structure_coords_before[:, 2]), -0.5 * params.pv_array.panel_thickness + np.amin(structure_coords_before[:, 2]), + params.pv_array.elevation - 0.5 * params.pv_array.panel_thickness, ) assert np.isclose( @@ -124,7 +129,8 @@ def __init__(self, domain, x_shift, y_shift, z_shift): np.amax(structure_coords_before[:, 1]), 0.5 * params.pv_array.panel_span ) assert np.isclose( - np.amax(structure_coords_before[:, 2]), 0.5 * params.pv_array.panel_thickness + np.amax(structure_coords_before[:, 2]), + params.pv_array.elevation + 0.5 * params.pv_array.panel_thickness, ) # Move the mesh by the amount prescribed in u_delta diff --git a/pvade/tests/test_solve.py b/pvade/tests/test_solve.py index 8734207d..68dfde89 100644 --- a/pvade/tests/test_solve.py +++ b/pvade/tests/test_solve.py @@ -37,6 +37,18 @@ def test_flow_3dpanels(): domain = FSIDomain(params) domain.read_mesh_files(rootdir + "/pvade/tests/input/mesh/panels3d/", params) + """ + for this mesh: + stream_rows: 7 + span_rows: 1 + elevation: 1.5 + stream_spacing: 7.0 + panel_chord: 2.0 + panel_span: 7.0 + panel_thickness: 0.1 + tracker_angle: -30.0 + """ + print("fluid shape = ", np.shape(domain.fluid.msh.geometry.x)) print("struct shape = ", np.shape(domain.structure.msh.geometry.x)) @@ -71,8 +83,12 @@ def test_flow_3dpanels(): print("max_pressure = ", max_pressure) assert not np.any(np.isnan(flow.p_k.x.array)) - max_velocity_truth = 18.205784057651652 - max_pressure_truth = 56.175743310367395 + # max_velocity_truth = 18.205784057651652 + # max_pressure_truth = 56.175743310367395 + + max_velocity_truth = 18.613277617512917 + max_pressure_truth = 73.50109146276495 + assert np.isclose(max_velocity, max_velocity_truth, rtol=rtol) assert np.isclose(max_pressure, max_pressure_truth, rtol=rtol) @@ -310,7 +326,7 @@ def test_fsi2(): # print('pos_data = ', pos_data) - assert np.allclose(pos_data, pos_data_truth) + assert np.allclose(pos_data[:, 0:2], pos_data_truth) print(lift_and_drag_data) # assert np.allclose(lift_and_drag_data[:, 0:3], lift_and_drag_data_truth[:, 0:3]) # needs new truth values to pass, mesh has changed diff --git a/pvade_main.py b/pvade_main.py index 8bc322aa..6b3d0462 100644 --- a/pvade_main.py +++ b/pvade_main.py @@ -16,6 +16,10 @@ import os from mpi4py import MPI +""" +to run: python -u pvade_main.py --input_file=input/duramat_case_study.yaml --domain.l_char=3.0 --pv_array.torque_tube_radius=0.1 --pv_array.torque_tube_separation=0.2 --pv_array.stream_rows=2 --pv_array.tracker_angle=[20,20] --pv_array.modules_per_span=10 +""" + def main(input_file=None): # Get the path to the input file from the command line @@ -38,6 +42,7 @@ def main(input_file=None): domain = FSIDomain(params) if params.general.input_mesh_dir is not None: domain.read_mesh_files(params.general.input_mesh_dir, params) + else: domain.build(params) @@ -60,7 +65,7 @@ def main(input_file=None): structure = Structure(domain, params) structure.build_boundary_conditions(domain, params) # # # Build the fluid forms - structure.build_forms(domain, params) + structure.build_forms(domain, params, flow) else: structure = None