""" =============================================================================== Diagnostic Script: PyFluent I/O and Physics Math Verification =============================================================================== Use this script to verify that Python can successfully: 1. Push 3D arrays to Fluent UDMs. 2. Trigger the Primal and Adjoint solvers. 3. Extract 3D gradient arrays from Fluent back into Python memory. =============================================================================== """ import os import numpy as np import scipy.ndimage as ndimage import ansys.fluent.core as pyfluent import matplotlib.pyplot as plt %matplotlib qt # ============================================================================= # 1. GLOBAL CONFIGURATION (Ensure these match your Fluent setup) # ============================================================================= NX, NY, NZ = 160, 80, 16 TOTAL_CELLS = NX * NY * NZ VOL_FRAC_TARGET = 0.50 K_FLUID, K_SOLID = 0.6, 205.0 ALPHA_MAX = 1e8 Q_ALPHA, Q_T = 0.1, 3.0 FILTER_RADIUS = 2.0 HEAVISIDE_BETA, HEAVISIDE_ETA = 1.0, 0.5 FLUID_ZONE_ID = 53 # IMPORTANT: Verify this ID in Fluent! CHIP_REPORT_DEF = "report-def-chip-avg-temp" # IMPORTANT: Verify exact spelling in Fluent! # ============================================================================= # 2. PHYSICS INTERPOLATION FUNCTIONS # ============================================================================= def extract_fluent_console_result(report_text): """ Extracts the last number (Net average) from a Fluent multi-line report string. """ # .split() automatically separates the string by any spaces, tabs, or newlines elements = report_text.split() # Grab the very last element [-1] and convert it to a float return float(elements[-1]) def string2double(file_name): #Converts the text (console-format) result files into double varibles. text=open(file_name,"r") double_number = text.readlines() double_number=double_number[-1] # Get last word in String double_number = double_number.split(' ')[-1] double_number = float(double_number) return double_number def calculate_brinkman_alpha(gamma): return ALPHA_MAX * ((gamma * (1.0 + Q_ALPHA)) / (gamma + Q_ALPHA)) def calculate_brinkman_derivative(gamma): return ALPHA_MAX * (Q_ALPHA * (1.0 + Q_ALPHA)) / ((gamma + Q_ALPHA)**2) def calculate_ramp_k(gamma): return K_FLUID + (gamma / (1.0 + Q_T * (1.0 - gamma))) * (K_SOLID - K_FLUID) def calculate_ramp_derivative(gamma): return (K_SOLID - K_FLUID) * ((1.0 + Q_T) / (1.0 + Q_T * (1.0 - gamma))**2) def apply_heaviside(gamma, beta, eta): num = np.tanh(beta * eta) + np.tanh(beta * (gamma - eta)) den = np.tanh(beta * eta) + np.tanh(beta * (1.0 - eta)) return num / den def compute_raw_sensitivities(gamma, v_p, v_adj, grad_T_p, grad_T_adj): v_dot_vadj = np.sum(v_p * v_adj, axis=1) sens_momentum = -1.0 * v_dot_vadj * calculate_brinkman_derivative(gamma) gradT_dot_gradTadj = np.sum(grad_T_p * grad_T_adj, axis=1) sens_thermal = -1.0 * gradT_dot_gradTadj * calculate_ramp_derivative(gamma) return sens_momentum + sens_thermal def generate_robust_pins(centroids, nx=160, ny=80, num_pins_x=20, num_pins_y=5, pin_width_cells=2): X = centroids[:, 0] Y = centroids[:, 1] # 1. Map physical X and Y exactly to integer grid indices (0 to 159, 0 to 79) x_min, x_max = X.min(), X.max() y_min, y_max = Y.min(), Y.max() idx_x = np.round((X - x_min) / (x_max - x_min) * (nx - 1)).astype(int) idx_y = np.round((Y - y_min) / (y_max - y_min) * (ny - 1)).astype(int) # 2. Create a clean 2D boolean array for the pattern is_solid_2d = np.zeros((nx, ny), dtype=bool) centers_x = np.linspace(0, nx - 1, num_pins_x, dtype=int) centers_y = np.linspace(0, ny - 1, num_pins_y, dtype=int) radius = pin_width_cells // 2 for cx in centers_x: for cy in centers_y: x_start = max(0, cx - radius) x_end = min(nx, cx + radius + (pin_width_cells % 2)) y_start = max(0, cy - radius) y_end = min(ny, cy + radius + (pin_width_cells % 2)) is_solid_2d[x_start:x_end, y_start:y_end] = True # 3. Map the clean 2D integer mask back onto the 3D unstructured cells natively is_solid_flat = is_solid_2d[idx_x, idx_y] return np.where(is_solid_flat, 1.0, 0.001) def generate_pins_from_centroids(centroids, num_pins_x=20, num_pins_y=5, pin_thickness=0.25): """ Evaluates a 2.5D sparse pin pattern directly on Fluent's unstructured cell coordinates. Args: centroids: NumPy array of shape (N, 3) containing X, Y, Z of every cell. num_pins_x, num_pins_y: Number of pins in each direction. pin_thickness: Fraction of the grid spacing the pin takes up (0.0 to 1.0). """ X = centroids[:, 0] Y = centroids[:, 1] # We ignore Z (centroids[:, 2]) entirely. This enforces the 2.5D extrusion constraint! # 1. Find the physical bounding box of the domain x_min, x_max = X.min(), X.max() y_min, y_max = Y.min(), Y.max() # 2. Calculate the pitch (spacing) for the pins pitch_x = (x_max - x_min) / num_pins_x pitch_y = (y_max - y_min) / num_pins_y # 3. Find the local coordinate within each pitch grid (normalized from 0.0 to 1.0) # Adding a tiny tolerance prevents edge-case wrapping issues local_x = ((X - x_min + 1e-9) % pitch_x) / pitch_x local_y = ((Y - y_min + 1e-9) % pitch_y) / pitch_y # 4. Define the solid pin boundaries (centered in the pitch cell) margin = (1.0 - pin_thickness) / 2.0 is_pin_x = (local_x >= margin) & (local_x <= (1.0 - margin)) is_pin_y = (local_y >= margin) & (local_y <= (1.0 - margin)) # 5. A cell is solid only if it falls inside both the X and Y pin boundaries is_solid = is_pin_x & is_pin_y # 6. Map to physical gamma values gamma_flat = np.where(is_solid, 1.0, 0.001) return gamma_flat # ============================================================================= # 3. SINGLE-PASS TEST FUNCTION # ============================================================================= def test_pyfluent_io(): print(">>> 1. Launching PyFluent...") solver = pyfluent.launch_fluent(show_gui=True, mode="solver", precision="double", processor_count=8) print(">>> 2. Reading Case File...") abs_case_path = r"C:\0-Work\0-Projects\260601_PM62Ext_Pinfin_CFD_adjoint_optim\04_scripting\260702_Fluent_setup_Topo_opt_pyfluent_v5.cas.h5" solver.file.read_case(file_name=abs_case_path) print(">>> 3. Generating Coordinate-Based Initial Pattern...") zone_name = "fluid_main" domain_name = "mixture" sv_data = solver.fields.solution_variable_data # --- 0. Initialize Fluent solver --- solver.settings.solution.initialization.initialization_type = "standard" solver.settings.solution.initialization.initialize() # --- 1. EXTRACT CENTROIDS --- centroids_dict = sv_data.get_data(variable_name="SV_CENTROID", zone_names=[zone_name], domain_name=domain_name) centroids = np.array(centroids_dict[zone_name]).reshape(-1, 3) num_cells = centroids.shape[0] # Exact ACTIVE cell count # --- 2. GENERATE PATTERN --- # Using the robust math for perfect sharp pins! #gamma_flat = generate_pins_from_centroids(centroids, num_pins_x=20, num_pins_y=7, pin_thickness=0.3) gamma_flat = generate_robust_pins(centroids, nx=NX, ny=NY, num_pins_x=5, num_pins_y=3, pin_width_cells=2) gamma_flat = apply_heaviside(gamma_flat, HEAVISIDE_BETA, HEAVISIDE_ETA) alpha_flat = calculate_brinkman_alpha(gamma_flat) k_flat = calculate_ramp_k(gamma_flat) solver.tui.define.user_defined.user_defined_memory(0) # Deallocate all memory solver.tui.define.user_defined.user_defined_memory(6) # Reallocate 6 fresh, zeroed blocks #solver.tui.define.user_defined.user_defined_memory(1) # Reallocate 6 fresh, zeroed blocks print(">>> 4. Pushing Data to Fluent UDMs (Pure 1D Native Stride)...") # --- 1. PULL THE NATIVE 1D SHAPE --- udm_dict = sv_data.get_data(variable_name="SV_UDM_I", zone_names=[zone_name], domain_name=domain_name) # Convert to a strict 1D numpy array arr_udm_flat = np.array(udm_dict[zone_name]).flatten() #arr_udm_flat = np.array(udm_dict[zone_name]) # --- 2. WIPE THE MEMORY --- # Create a fresh 1D array of the exact same size (wipes all 1e8 ghosts) arr_udm_flat = np.zeros_like(arr_udm_flat) # --- 3. 1D STRIDE ASSIGNMENT (Interleaved Mapping) --- # arr_udm_flat[0::6] accesses all UDM-0 slots. # We slice [:num_cells] to ensure we only fill the active cells and leave ghosts as 0.0 arr_udm_flat[0::6][:num_cells] = gamma_flat # UDM-0 arr_udm_flat[1::6][:num_cells] = alpha_flat # UDM-1 arr_udm_flat[2::6][:num_cells] = k_flat # UDM-2 # --- 4. PUSH DATA --- # Pass the pure 1D NumPy array directly (NO .flatten(), NO .tolist()) sv_data.set_data( variable_name="SV_UDM_I", zone_names_to_data={zone_name: arr_udm_flat}, #zone_names_to_data={zone_name: gamma_flat}, domain_name=domain_name ) solver.tui.display.re_render() # Force Fluent to keep temporary solver memory (gradients) alive solver.tui.solve.set.expert("yes", "yes", "yes", "yes") try: solver.settings.solution.controls.advanced.expert = {'keep_temporary_memory' : True} except: pass print("Solving Primal Equations...") solver.settings.solution.run_calculation.iterate(iter_count=10) print(">>> 6. Evaluating Report Definition...") try: t_avg_list = solver.settings.solution.report_definitions.compute(report_defs=[CHIP_REPORT_DEF]) # Robust dictionary extraction depending on PyFluent version if isinstance(t_avg_list, list): t_avg_celsius = t_avg_list[0][CHIP_REPORT_DEF][0] - 273.15 else: t_avg_celsius = t_avg_list[CHIP_REPORT_DEF][0] - 273.15 print(f" [SUCCESS] Extracted Chip Temp: {t_avg_celsius:.2f} C") except Exception as e: print(f" [FAILED] Could not read report def '{CHIP_REPORT_DEF}'. Error: {e}") print(">>> 7. Running Adjoint Solver...") try: solver.tui.adjoint.run.initialize() solver.tui.adjoint.run.initialize_strategy() solver.tui.adjoint.run.iterate(15) print(" [SUCCESS] Adjoint solve completed.") except Exception as e: print(f" [FAILED] Adjoint solver failed. Is it set up in the GUI? Error: {e}") print(">>> 8. Extracting Fields via Solution Variables...") sv_data = solver.fields.solution_variable_data def get_sv_np(sv_name_str): data_dict = sv_data.get_data(variable_name=sv_name_str, zone_names=[zone_name], domain_name=domain_name) return np.array(data_dict[zone_name]).flatten()[:num_cells] try: # 1. Primal Velocity v_p = np.column_stack((get_sv_np('SV_U'), get_sv_np('SV_V'), get_sv_np('SV_W'))) # 2. Primal Temperature Gradients (Read from custom UDMs!) udm_dict = sv_data.get_data(variable_name='SV_UDM_I', zone_names=[zone_name], domain_name=domain_name) udm_master_flat = np.array(udm_dict[zone_name]).flatten() # Reshape to SEQUENTIAL format based on total allocated memory num_allocated_ext = udm_master_flat.size // 6 udm_master_2d = udm_master_flat.reshape(6, num_allocated_ext) # EXTRACT ROWS 3, 4, AND 5 (and drop the ghost cells!) dT_x = udm_master_2d[3, :num_cells] dT_y = udm_master_2d[4, :num_cells] dT_z = udm_master_2d[5, :num_cells] grad_T_p = np.column_stack((dT_x, dT_y, dT_z)) # 3. Adjoint Velocity v_adj = np.column_stack((get_sv_np('SV_ADJOINT_U'), get_sv_np('SV_ADJOINT_V'), get_sv_np('SV_ADJOINT_W'))) # 4. Adjoint Temperature Gradients adj_t_g_dict = sv_data.get_data(variable_name='SV_ADJOINT_T_G', zone_names=[zone_name], domain_name=domain_name) adj_t_g_flat = np.array(adj_t_g_dict[zone_name]).flatten() grad_T_adj = adj_t_g_flat.reshape(-1, 3)[:num_cells, :] print(f" [SUCCESS] Fields extracted! v_p shape: {v_p.shape} matches gamma: {gamma_flat.shape}") except Exception as e: print(f" [FAILED] Extraction error: {e}") print(">>> 9. Testing Sensitivity Math...") try: sens_3d_flat = compute_raw_sensitivities(gamma_flat, v_p, v_adj, grad_T_p, grad_T_adj) sens_3d = sens_3d_flat.reshape((NX, NY, NZ)) sens_2d = np.mean(sens_3d, axis=2) sens_2d_filtered = ndimage.gaussian_filter(sens_2d, sigma=FILTER_RADIUS) print(" [SUCCESS] Sensitivity arrays processed and filtered.") except Exception as e: print(f" [FAILED] Array math failed. Check mesh dimensions NX/NY/NZ. Error: {e}") print(">>> 10. Test Complete. Closing Fluent.") solver.exit() if __name__ == "__main__": test_pyfluent_io() #%% #def visualize_gamma_3d(centroids, gamma_flat): """ Creates an interactive 3D scatter plot of the solid cells. """ print(">>> Generating 3D Interactive Plot...") # Filter for only the solid cells (pins) solid_mask = gamma_flat > 0.5 solid_coords = centroids[solid_mask] if len(solid_coords) == 0: print(" [WARNING] No solid cells found! The pattern generation failed.") #return print(f" [INFO] Plotting {len(solid_coords)} solid cells out of {len(centroids)} total.") # Create the 3D plot fig = plt.figure(figsize=(10, 8)) ax = fig.add_subplot(111, projection='3d') # Plot the centroids of the solid cells # s=1 is the marker size, alpha=0.5 adds slight transparency ax.scatter(solid_coords[:, 0], solid_coords[:, 1], solid_coords[:, 2], c='red', marker='s', s=2, alpha=0.8) # Formatting ax.set_title("Python-Side Gamma Initialization (Solid Cells Only)") ax.set_xlabel("X Coordinate") ax.set_ylabel("Y Coordinate") ax.set_zlabel("Z Coordinate") # Ensure axes scales are equal so the geometry doesn't look stretched # Note: set_aspect('equal') in 3D can be finicky depending on matplotlib version, # so we set the limits manually based on the domain bounding box max_range = np.array([ centroids[:,0].max()-centroids[:,0].min(), centroids[:,1].max()-centroids[:,1].min(), centroids[:,2].max()-centroids[:,2].min() ]).max() / 2.0 mid_x = (centroids[:,0].max()+centroids[:,0].min()) * 0.5 mid_y = (centroids[:,1].max()+centroids[:,1].min()) * 0.5 mid_z = (centroids[:,2].max()+centroids[:,2].min()) * 0.5 ax.set_xlim(mid_x - max_range, mid_x + max_range) ax.set_ylim(mid_y - max_range, mid_y + max_range) ax.set_zlim(mid_z - max_range, mid_z + max_range) # Launch the interactive window plt.show()