From cdfa6c51584d1b699373a57b3cd79bbcb01b153d Mon Sep 17 00:00:00 2001 From: Bryce Mazurowski Date: Wed, 16 Apr 2025 10:54:45 -0500 Subject: [PATCH 1/7] fix: update skfem version spec to avoid issue with default mutable datatype --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index d599151..731d36b 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,5 @@ numpy -scipy==1.9.2 +scipy>=1.9.2 h5py vtk pytest From e6ad15462547934a795e0f4921f075b954b1d5b6 Mon Sep 17 00:00:00 2001 From: Bryce Mazurowski Date: Wed, 16 Apr 2025 10:57:12 -0500 Subject: [PATCH 2/7] feat: Added initial guess to flowpath solver Save as member of ThermohydraulicsThermalSolver after successful solve Then pass as arg to flowpath.solve() --- srlife/thermal.py | 5 ++++- srlife/thermohydraulics/flowpath.py | 12 ++++++++++-- 2 files changed, 14 insertions(+), 3 deletions(-) diff --git a/srlife/thermal.py b/srlife/thermal.py index 14f21c6..3c4243a 100644 --- a/srlife/thermal.py +++ b/srlife/thermal.py @@ -86,6 +86,7 @@ def __init__(self, pset=solverparams.ParameterSet()): self.solid_params = pset.get_default("solid", solverparams.ParameterSet()) self.thermo_params = pset.get_default("fluid", solverparams.ParameterSet()) + self.conv_fluid_nodal_T = np.array([]) def solve_receiver( self, @@ -311,7 +312,9 @@ def solve_fluid(self, i, time, dt): ) # Solve for the fluid temperatures - nodal_temps = model.solve(time) + nodal_temps = model.solve(time, self.conv_fluid_nodal_T) + # save these values to use as initial guess in next solve + self.conv_fluid_nodal_T = nodal_temps # Update the stored fluid temperature and flow velocities flow_rates, tube_temperatures = model.recover_tube_results( diff --git a/srlife/thermohydraulics/flowpath.py b/srlife/thermohydraulics/flowpath.py index 40a17b0..1b11bf0 100644 --- a/srlife/thermohydraulics/flowpath.py +++ b/srlife/thermohydraulics/flowpath.py @@ -275,6 +275,11 @@ def __init__( self.miter = 100 self.verbose = verbose + # Save previous solution for use as initial guess + # to start, this is empty + self.T_prev_sol = np.array([]) + self.has_prev_conv_T = False + def add_panel(self, weights, ri, h, metal_temp, material): """ Construct and add the standard panel -> manifold link @@ -342,7 +347,7 @@ def _setup(self): self.dof_map.append(list(range(self.nvals, self.nvals + obj.size))) self.nvals += obj.size - def solve(self, t): + def solve(self, t, T_init_guess=np.array([])): """ Solve for the current fluid temperatures @@ -353,7 +358,10 @@ def solve(self, t): self.t = t # Significant decision... - T = np.zeros((self.nvals)) + if T_init_guess.size != 0: + T = T_init_guess + else: + T = np.zeros((self.nvals)) # Initial residual R, J = self.RJ(T) From faa7d7af6d354582c4415d6088637c01bf441b80 Mon Sep 17 00:00:00 2001 From: Bryce Mazurowski Date: Wed, 16 Apr 2025 10:58:29 -0500 Subject: [PATCH 3/7] fix: updated deprecated np.trapz to np.trapezoid --- srlife/damage.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/srlife/damage.py b/srlife/damage.py index 4b9053e..56043e5 100644 --- a/srlife/damage.py +++ b/srlife/damage.py @@ -1033,7 +1033,7 @@ def calculate_volume_flaw_flattened_eq_stress( else: sigma_e[sigma_e < 0] = 0 g = ( - np.trapz( + np.trapezoid( (sigma_e / (sigma_e_max + self.tolerance)) ** Navg[..., None, None], time, axis=0, @@ -1196,7 +1196,7 @@ def calculate_surface_flaw_flattened_eq_stress( else: sigma_e[sigma_e < 0] = 0 g = ( - np.trapz( + np.trapezoid( (sigma_e / (sigma_e_max + self.tolerance)) ** Navg[..., None, None, None], time, @@ -1474,7 +1474,7 @@ def calculate_volume_flaw_element_log_reliability( # Calculating ratio of cyclic stress to max cyclic stress (in one cycle) # for time-independent and time-dependent cases g = ( - np.trapz( + np.trapezoid( (pstress / (pstress_max + 1.0e-14)) ** Navg[..., None], time, axis=0, @@ -1556,7 +1556,7 @@ def calculate_surface_time_dep_pstress( # Calculating ratio of cyclic stress to max cyclic stress (in one cycle) # for time-independent and time-dependent cases g = ( - np.trapz( + np.trapezoid( (surf_pstress / (surf_pstress_max + 1.0e-14)) ** Navg[..., None, None], time, axis=0, @@ -1743,7 +1743,7 @@ def calculate_volume_flaw_avg_normal_stress( sigma_n_0 = sigma_n else: g = ( - np.trapz( + np.trapezoid( (sigma_n / (sigma_n_max + self.tolerance)) ** Navg[..., None, None], time, axis=0, @@ -1890,7 +1890,7 @@ def calculate_surface_flaw_time_dep_normal_stress( sigma_n_0 = sigma_n else: g = ( - np.trapz( + np.trapezoid( (sigma_n / (sigma_n_max + self.tolerance)) ** Navg[..., None, None, None], time, From de1e2e6b8a2cd996fd070a6dd4b4267368b6aa95 Mon Sep 17 00:00:00 2001 From: Bryce Mazurowski Date: Wed, 16 Apr 2025 11:39:59 -0500 Subject: [PATCH 4/7] feat: Added script interface to srlife to allow full analysis srlife/interface.py has funcitons to facilitate interacting with library Added suite of tests to check funcitonality Added manual workflows to test Updated test to work on any machine More ssmall fixes LastFix --- .github/workflows/formatting.yml | 2 +- .github/workflows/test.yml | 2 +- srlife/interface.py | 1023 ++++++++++++++++++++++++++++++ srlife/receiver.py | 2 +- test/Rec_Gold.hdf5 | Bin 0 -> 115736 bytes test/data_6_d20_hr6.csv | 32 + test/data_6_d20_hr7.csv | 32 + test/test_interface.py | 551 ++++++++++++++++ 8 files changed, 1641 insertions(+), 3 deletions(-) create mode 100644 srlife/interface.py create mode 100644 test/Rec_Gold.hdf5 create mode 100644 test/data_6_d20_hr6.csv create mode 100644 test/data_6_d20_hr7.csv create mode 100644 test/test_interface.py diff --git a/.github/workflows/formatting.yml b/.github/workflows/formatting.yml index 480c3bf..2f02d52 100644 --- a/.github/workflows/formatting.yml +++ b/.github/workflows/formatting.yml @@ -1,5 +1,5 @@ name: Check code formatting -on: [push, pull_request] +on: [push, pull_request, workflow_dispatch] jobs: check-formatting: runs-on: ubuntu-latest diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 6b17bca..6ec1043 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -1,5 +1,5 @@ name: tests -on: [push, pull_request] +on: [push, pull_request, workflow_dispatch] jobs: test-linux: runs-on: ubuntu-latest diff --git a/srlife/interface.py b/srlife/interface.py new file mode 100644 index 0000000..0a02b4c --- /dev/null +++ b/srlife/interface.py @@ -0,0 +1,1023 @@ +""" + This module has all of the helper functions needed to run an srlife analysis. + It is meant to allow a simple script interface to define a problem and run the + analysis. + + One such example is defined in test/test_interface.py +""" + +import numpy as np +from scipy.interpolate import RegularGridInterpolator + +# srlife specific functions +from srlife import ( + receiver, + solverparams, +) + + +# UNIT CONVERSION FUNCTIONS +def convert_mm_to_m(mm_in): + """ + Unit conversion from mm to m + + Args: mm_in (double): qty to convert + """ + return mm_in / 1000 + + +def convert_m_to_mm(m_in): + """ + Unit conversion from m to mm + + Args: m_in (double): qty to convert + """ + return m_in * 1000 + + +def convert_kWm2_to_Wmm2(kWm2_in): + """ + Unit conversion from kW/m^2 to W/mm^2 + + Args: kWm2_in (double): qty to convert + """ + return kWm2_in * 1000e-6 + + +def convert_Pa_to_MPa(Pa_in): + """ + Unit conversion from Pa to MPa + + Args: Pa_in (double): qty to convert + """ + return Pa_in * 1e-6 + + +def convert_MPa_to_Pa(MPa_in): + """ + Unit conversion from MPa to Pa + + Args: MPa_in (double): qty to convert + """ + return MPa_in * 1e6 + + +def convert_C_to_K(C_in): + """ + Unit conversion from C to K + + Args: C_in (double): qty to convert + """ + return C_in + 273.15 + + +def read_month_day_hour_flux_file(flux_file_name, data_shape, z_offset_in): + """ + Read and process the flux data from a solarPILOT file + + Args: + flux_file_name (string): path and name of flux file from solarpilot + data_shape (list[int]): list with n_rows and n_cols of flux data + z_offset_in (double): z offset between srlife model and solarpilot model + this will shift z data that is read up or down + + Returns: + z_data (np.array(double)): vector with z positions from flux sampling + theta_data (np.array(double)): vector with theta positions from flux sampling + flux_data (np.array(double)): matrix with flux data at sampled positions + """ + # read flux data + # massage into format + # output + row, col = data_shape + z_data = np.zeros([row, 1]) + theta_data = np.zeros([1, col]) + flux_data = np.zeros([row, col]) + try: + theta_data = np.loadtxt( + flux_file_name, delimiter=",", skiprows=6, usecols=range(1, col + 1) + )[0] + z_data = np.loadtxt( + flux_file_name, delimiter=",", skiprows=7, usecols=range(0, 1) + ) + z_data += z_offset_in + flux_data = np.loadtxt( + flux_file_name, delimiter=",", skiprows=7, usecols=range(1, col + 1) + ) + # Data is on cylinder, so 180 will repeat + theta_data = np.append(theta_data, -180) + flux_data = np.hstack((flux_data, flux_data[:, 0][:, None])) + except FileNotFoundError: + # if we cannot find the file + # print a message and just return all zeros + print( + f"File could not be found!!! Returning all zeros\n filename: {flux_file_name}" + ) + # convert height to mm + z_data = convert_m_to_mm(z_data) + return z_data, theta_data, flux_data + + +def get_flux_interpolators_from_data_files( + times, start_time, end_time, month, day, flux_data_dir, flux_data_shape, z_offset +): + """ + Read flux files and create interpolator functions for each month, day, hour + + Args: + times (list(int)): list of index of hours used for analysis + Could probably remove this if (end_time - start_time) == len(times) + start_time (int): first hour to read flux data + end_time (int): last hour to read flux data + month (int): month to read flux data + day (int): day to read flux data + flux_data_dir (string): path to flux data files + flux_data_shape (list(int)): n_rows and n_cols of flux data in flux files + z_offset_in (double): z offset between srlife model and solarpilot model + this will shift z data that is read up or down + + Returns: + flux_interpolators_by_hour (list(RegularGridInterpolators)): list with the interpolating + functions for each hour in times. + """ + print("Making flux interpolators from data") + # Time zero is NOT analyzed, so we do not make a fluxFn for that + flux_interpolators_by_hour = [0] * (len(times) - 1) + for i_count, i_time in enumerate(range(start_time, end_time + 1)): + flux_filename = f"{flux_data_dir}/data_{month}_d{day}_hr{int(i_time)}.csv" + z_pts, theta_pts, flux_pts = read_month_day_hour_flux_file( + flux_filename, flux_data_shape, z_offset + ) + interpFun = RegularGridInterpolator( + (z_pts, theta_pts), + flux_pts, + method="linear", + bounds_error=False, + fill_value=0.0, + ) + flux_interpolators_by_hour[i_count] = interpFun + return flux_interpolators_by_hour + + +def apply_tube_flux_bcs( + tube, + tube_absorbance, + tube_rec_theta, + flux_interpolators_by_hour, +): + """ + Actually apply the heat flux to the tube across width and length + + Args: + tube (srlife.Tube): srlife object for this tube + tube_z (list[double]): list of z coordinates of tube nodes + tube_theta (list[double]): list of theta coordinates of tube nodes + tube_absorbance (double): absorbance of tubes + tube_rec_theta (double): receiver theta value for this tube + flux_interpolators_by_hour (list[RegularGridInterpolator]): + list of interpolating functions for heat flux on tubes + """ + num_steps = len(flux_interpolators_by_hour) + tube_theta = np.linspace(0, 2 * np.pi, tube.nt + 1)[:-1] + tube_z = np.linspace(0, tube.h, tube.nz) + # add extra time step for time zero + tube_flux = np.zeros([num_steps + 1, len(tube_theta), len(tube_z)]) + for i_hour in range(num_steps): + # for each hour of flux data + # sample tube at tube_rec_theta and z + for i_z, z in enumerate(tube_z): + # for each zPt + pt_z_theta = [z, tube_rec_theta] + rec_flux = convert_kWm2_to_Wmm2( + flux_interpolators_by_hour[i_hour](pt_z_theta) + ) + # apply across crown so only sunSide gets flux + for i_theta, theta in enumerate(tube_theta): + # for each thetaPt around tube + if theta > 270 or theta < 90: + tube_flux[i_hour, i_theta, i_z] = rec_flux[0] * np.cos( + theta * np.pi / 180 + ) + else: + tube_flux[i_hour, i_theta, i_z] = 0.0 + # make bc object for tube with heat flux data + # NOTE: this is done here and not in cerate_tube because if we initialize + # bc data to np.zeros, flux will always be zero unless we reload from a file + tube.set_bc(apply_heatflux_bc(tube, tube_flux * tube_absorbance, "outer"), "outer") + # while doing this, also apply convective BCs to inner face + tube.set_bc(apply_convective_bc(tube, "inner"), "inner") + + +def apply_convective_bc(tube, loc): + """ + Defines and applies a convective BC to a surface of a tube + specificed by loc + + Args: + tube (srlife.Tube): srlife object for this tube + loc (string): "inner" or "outer" to specify which surface BC is applied to + + """ + nz = tube.nz + if loc == "outer": + r = tube.r + else: + r = tube.r - tube.t + return receiver.ConvectiveBC( + r, tube.h, nz, tube.times, np.zeros((len(tube.times), nz)) + ) + + +def apply_heatflux_bc(tube, data, loc): + """ + Defines and applied a heat flux BC to a surface of a tube object + specified by loc + + Args: + tube (srlife.Tube): srlife object for this tube + data (list[double]): 3d list of heat flux bcs [hour, theta, z] + loc (string): "inner" or "outer" to specify which surface BC is applied to + + """ + if loc == "outer": + r = tube.r + else: + r = tube.r - tube.t + nt = tube.nt + nz = tube.nz + return receiver.HeatFluxBC(r, tube.h, nt, nz, tube.times, data) + + +def calc_and_write_tube_flux_bcs( + rec, + ass_tube_per_panel, + num_panels, + flux_interpolators_by_hour, + tube_absorbance, +): + """ + Applies flux boundary conditions to analyzed tubes. Flux value is interpolated based on + centerline tube location. It is applied along the tube face as a cosine distribution + with max at sun-facing centerline to zero at 90 deg in either direction. + + Args: + rec (Receiver): the receiver object that the tubes belong to + ass_tube_per_panel (int): the assumed number of tubes per panel, or analyzed number of tubes + num_panels (int): number of panels in the receiver #BPMToDo: remove this and get from rec + flux_interpolators_by_hour (list(RegularGridInterpolators)): list of flux interpolators + on receiver by hour of analysis + tube_z (list(double)): z nodes points along tube + tube_theta (list(double)): theta nodes points along tube + tube_absorbance (double): absorbance of tubes + write_tube_flux_to_csv (bool): flag to dump flux data per tube to a csv file to debug + + Returns: None + """ + print("Assigning flux BCs to tubes") + for panel_key, panel in rec.panels.items(): + # in each panel of rec + for tube_key, tube in panel.tubes.items(): + print(f"Panel: {panel_key} Tube: {tube_key}") + # in each tube of panel + tube_rec_theta = 180 - (360 * int(panel_key) / num_panels) + if ass_tube_per_panel == 1: + # if only one tube, place at center of panel + tube_rec_theta += 360 / (2 * num_panels) + else: + # if more than one tube, evenly space them including edges + tube_rec_theta -= ( + 360 / num_panels * int(tube_key) / (ass_tube_per_panel - 1) + ) + if tube_rec_theta == -180: + # wrap data around tube + tube_rec_theta = 180 + apply_tube_flux_bcs( + tube, + tube_absorbance, + tube_rec_theta, + flux_interpolators_by_hour, + ) + + +def set_rec_flow_paths(rec, panel_flow_path, mass_flow_per_path, T_in_per_path): + """ + Set flowpaths within a receiver object + + Args: + rec (Receiver): the receiver object to set flowpaths for + panel_flow_path (list[string]): panel names indicating the flowpath. Can have a list + of lists if more that one flowpath + mass_flow_per_path (list[double]): list with entries of initial mass flow rate for + each path + T_in_per_path (list[double]): list with entries of inlet temp for each path + """ + print("Setting receiver flowpaths") + T_in_per_path = convert_C_to_K(T_in_per_path) + # add flowpaths to receiver + for i_path, flow_path_panels in enumerate(panel_flow_path): + # for each path in given list + print(f"Path {i_path}") + flow_path = [0] * len(flow_path_panels) + for i_panel, panel in enumerate(flow_path_panels): + # for each panel in path + print(f"Panel: {panel}") + flow_path[i_panel] = panel + panel = rec.panels[panel] + for tube in panel.tubes.values(): + # for each tube in panel + times = tube.times + data_shape = tube.outer_bc.data[:, 0, 0].shape + # initialize mass flow for path + mass_flow = mass_flow_per_path[i_path] + # convert mass flow from kg/s to kg/hr + mass_flow *= 3600 + mass_flow *= np.ones(data_shape) + # set inlet temp + T_in_flow_path = T_in_per_path[i_path] * np.ones_like(times) + rec.add_flowpath(flow_path, times, mass_flow, T_in_flow_path) + + +def check_outlet_temps_and_step_mass_flow( + rec, path, path_key, T_out_target, pct_err_outlet_temp, breaker +): + """ + Check outlet temp of all tubes from last panel in flowpath. + Average their value and use it to step mass flow rate. + + Args: + rec (srlife.Receiver): Receiver object who owns the flowpath we are + optimizing + path (OrderedDict): dictionary of flowpath information + path_key (string): name of flowpath + T_out_target (double): target outlet temperature + pct_err_outler_temp (double): acceptable error % for outlet temp (for convergence) + breaker (list[bool]): check whether all flowpaths have converged + """ + # All tubes in last panel go to manifold, so estimate + # temp as avg of all tubes output + panels = path.get("panels") + T_in_path = path["inlet_temp"][1:] + last_panel = rec.panels[panels[-1]] + num_tubes = last_panel.ntubes + num_time_steps = next(iter(last_panel.tubes.values())).ntime + T_out_path = np.zeros([num_tubes, num_time_steps - 1]) + for i_tube, tube in enumerate(last_panel.tubes.values()): + # for tubes in the last panel of the flow path + # results for tube are (time, nz) shape + # want all time except 0 and the last zPos + T_out_path[i_tube] = tube.axial_results.get("fluid_temperature")[1:, -1] + T_out_path = T_out_path.mean(axis=0) + print(f"T_out = {T_out_path - 273.15}") + if ( + any(T_out_path > T_out_target * (1 + pct_err_outlet_temp / 100)) is False + and any(T_out_path < T_out_target * (1 - pct_err_outlet_temp / 100)) is False + ): + # we have converged this flowpath + breaker[int(path_key)] = True + else: + mass_flow = path.get("mass_flow") + mass_flow[1:] = ( + mass_flow[1:] * (T_out_path - T_in_path) / (T_out_target - T_in_path) + ) + mass_flow[0] = mass_flow[1] + rec.flowpaths[path_key]["mass_flow"] = mass_flow + print(f"Mass flow: {mass_flow/3600}") + + +def optimize_mass_flow_rate_per_path( + rec, rec_filename, T_out_target, pct_err_outlet_temp, solver, save_heat_to_vtu +): + """ + Iteratively calculate the ideal mass flow rate for each path such that the + outlet temperature matches a target temperature. + + Args: + rec (Receiver): srlife Receiver object we are optimizing flow for + rec_filename (string): file name for the receiver object to optimize + T_out_target (double): target output temperature + pct_err_outlet_temp (double): acceptable percentage of error for the outlet temp + solver (managers.SolutionManager): srlife solver object + save_heat_to_vtu (bool): whether or not to save vtu files of results + """ + print("Optimizing mass flow rate to match outlet temp") + T_out_target = convert_C_to_K(T_out_target) + # make sure receiver is up to date + solver.receiver = rec + # Mass flow optimization + # limit of error on outlet temp + max_mass_flow_iter = 5 + for i_opt in range(max_mass_flow_iter): + print(f"Iteration = {i_opt}") + breaker = [False] * len(rec.flowpaths) + if i_opt > 0: + for path_key, path in rec.flowpaths.items(): + print(f"FlowPath: {path_key}") + # check last panel tube outlet temps + check_outlet_temps_and_step_mass_flow( + rec, path, path_key, T_out_target, pct_err_outlet_temp, breaker + ) + + if all(breaker): + print(f"Converged!!! solved in {i_opt} iterations") + if save_heat_to_vtu: + for panel_key, panel in rec.panels.items(): + for tube_key, tube in panel.tubes.items(): + tube.write_vtk(f"ht-tube-{panel_key}-{tube_key}") + break + # solve heat transfer problem + solver.solve_heat_transfer() + filename = rec_filename + "_ht_iter_" + str(i_opt) + ".hdf5" + rec.save(filename) + + +def calc_fluid_velocity(tube_mass_flow, fluid_rho, tube_Dh): + """ + Calculate the velocity of the fluid in a tube based on mass flow rate + + Args: + tube_mass_flow (double): mass flow rate in tube + fluid_rho (double): fluid density + tube_Dh (double): hydraulic diameter of tube + """ + return tube_mass_flow / fluid_rho * (4.0 / (np.pi * tube_Dh**2)) + + +def calc_reynolds_number(fluid_rho, fluid_vel, fluid_mu, tube_Dh): + """ + Calculate the Reynolds number for tube flow + + Args: + fluid_rho (double): fluid density + fluid_vel (double): fluid velocity + fluid_mu (double): fluid dynamic viscosity + tube_Dh (double): tube hydraulic diameter + + """ + return fluid_rho * fluid_vel * tube_Dh / fluid_mu + + +def calc_fd(Re, tube_eta, tube_Dh): + """ + Calculate Darcy friction factor from: + Zigrange and Sylvester 1985, A review of explicit friction factor equations. + J of Energy Resources Technology) + Equation 13 + + Args: + Re (double): Reynolds number + tube_eta (double): tube roughness + tube_Dh (double): tube hydraulic diameter + Returns: + fd (double): darcy friciton factor + """ + if Re < 4000: + # laminar flow + fd = 64 / Re + else: + fd = ( + 1 + / ( + -2 + * np.log10( + (tube_eta / tube_Dh) / 3.7 + - 5.02 / Re * np.log10((tube_eta / tube_Dh) / 3.7 + 13 / Re) + ) + ) + ** 2 + ) + return fd + + +def calc_friction_p_loss(fluid_rho, fluid_mu, tube_mass_flow, tube_dict, delta_l): + """ + Calculate pressure loss from friction in tubes + + Args: + fluid_rho (double): fluid density in tubes + fluid_mu (double): fluid dynamic visco in tubes + tube_mass_flow (double): mass flow rate in tube + tube_dict (dict): dictionary of tube specs + delta_l (double): length of tube + + Returns: + friction_p_loss (double): pressure loss from friction + """ + # tube roughness + tube_eta = tube_dict["eta"] + tube_Dh = convert_mm_to_m(tube_dict["od"] - 2 * tube_dict["t"]) + fluid_vel = calc_fluid_velocity(tube_mass_flow, fluid_rho, tube_Dh) + Re = calc_reynolds_number(fluid_rho, fluid_vel, fluid_mu, tube_Dh) + fd = calc_fd(Re, tube_eta, tube_Dh) + return fd * (0.5 * fluid_rho) * (fluid_vel**2) / tube_Dh * delta_l + + +def calc_manifold_bend_loss( + max_mass_flow, manifold_dict, manifold_rho, manifold_mu, num_bends +): + """ + Calculate pressure loss from flow through bends in manifold tubes + + Args: + max_mass_flow (double): mass flow in manifold + manifold_dict (dict): dictionary with manifold tube specs + manifold_rho (double): fluid density in manifold + manifold_mu (double): dynamic viscosity in manifold + num_bends (int): total number of bends in manifold tubes + """ + manifold_eta = manifold_dict["eta"] + manifold_pipe_t = convert_mm_to_m(manifold_dict["t"]) + manifold_pipe_id = convert_mm_to_m(manifold_dict["od"] - 2 * manifold_pipe_t) + manifold_pipe_bend_radius = convert_mm_to_m(manifold_dict["bend_radius"]) + manifold_eta = manifold_dict["eta"] + manifold_Dh = manifold_pipe_id + if manifold_eta / manifold_Dh > 0.0001: + K0 = 0.42 + elif manifold_eta / manifold_Dh == 0: + K0 = 0.21 + else: + K0 = 0.21 * (1 + 1000 * manifold_eta / manifold_Dh) + + manifold_vel = calc_fluid_velocity(max_mass_flow, manifold_rho, manifold_Dh) + Re = calc_reynolds_number(manifold_rho, manifold_vel, manifold_mu, manifold_Dh) + manifold_fd = calc_fd(Re, manifold_eta, manifold_Dh) + K = ( + K0 ** ((manifold_pipe_id / manifold_pipe_bend_radius) ** 0.5) + + 0.5 * np.pi * (manifold_pipe_bend_radius / manifold_pipe_id) * manifold_fd + ) + return num_bends * K * manifold_rho * manifold_vel**2 + + +def calc_path_total_p_loss( + rec, + path, + path_key, + delta_T, + tube_dict, + manifold_dict, + fluid, + rec_circumference, + num_panels, +): + """ + Calculate the total pressure loss in a flowpath from flow friction, + head loss, and flow through pipe bends. + + Args: + rec (Receiver): receiver object to analyze + path (dict): flowpath to get pressure loss through + path_key (string): name of flowpath + delta_T (double): temp step through flowpath + tube_dict (dict): dictionary of usefule tube info + manifold_dict (dict): dict of useful info for manifold tubes + fluid (thermalfluid.ThermalFluidMaterial): working fluid in receiver + rec_circumference (double): circumference of receiver + num_panels (int): number of panels in receiver + + Returns: + total_p_loss (double): total pressure loss in flowpath + """ + g = 9.81 # m/s^2 + # NOTE: this function works mostly with meters + act_tube_per_panel = tube_dict["ass_tube_per_panel"] * tube_dict["tube_mult"] + num_bends = num_panels * manifold_dict["bends_per_panel"] + + inlet_temp = path["inlet_temp"][0] + # get outlet avg outlet temp from each tube in last panel + # then average down again to single value + outlet_temp = np.mean( + [ + tube.axial_results.get("fluid_temperature")[:, -1][1:-1].mean() + for tube_key, tube in rec.panels[path["panels"][-1]].tubes.items() + ] + ) + max_mass_flow = path["mass_flow"].max() / 3600 # kg/s + tube_mass_flow = max_mass_flow / act_tube_per_panel + flow_path_length = convert_mm_to_m(tube_dict["h"]) * len(path["panels"]) + temp_steps = np.arange(inlet_temp, outlet_temp, delta_T) + rho_with_T = fluid.rho(temp_steps) * 1e9 # kg/m^3 + mu_with_T = fluid.mu(temp_steps) * 1000 / 3600 # Pa-s + + # Friction Loss + # pos spacing of temp change on tubes + delta_l = flow_path_length / np.size(temp_steps) + friction_p_loss = sum( + calc_friction_p_loss(rho, mu, tube_mass_flow, tube_dict, delta_l) + for rho, mu in zip(rho_with_T, mu_with_T) + ) + # Head Loss + avg_temp = 0.5 * (inlet_temp + outlet_temp) + head_p_loss = fluid.rho(avg_temp) * g * convert_mm_to_m(tube_dict["h"]) + # Manifold Loss + manifold_friction_p_loss = calc_friction_p_loss( + fluid.rho(avg_temp) * 1e9, + fluid.mu(avg_temp) * 1000 / 3600, + max_mass_flow, + manifold_dict, + (len(path["panels"]) - 1) * (rec_circumference / num_panels), + ) + manifold_bend_loss = calc_manifold_bend_loss( + max_mass_flow, + manifold_dict, + fluid.rho(avg_temp) * 1e9, + fluid.mu(avg_temp) * 1000 / 3600, + num_bends, + ) + return friction_p_loss + head_p_loss + manifold_bend_loss + manifold_friction_p_loss + + +def calc_p_loss_from_flows_temps(rec, tube_dict, manifold_dict, fluid, outlet_p): + """ + Calculate the pressure loss across the receiver based on the results + of a flow simulation. This captures friction loss, head loss, and bend losses + in the tubes and manifold pipes of the receiver + + Args: + rec (Receiver): the receiver object to calculate pressure loss + tube_dict (dict): a dictionary containing info about rec tubes + manifold_dict (dict): a dictionary containing info about manifold tubes + fluid (thermalfluid.ThermalFluidMaterial): thermal material of working fluid + outlet_p (double): pressure at flowpath outlets + """ + print("Calculating pressure loss from flow at temps") + # receiver details + num_panels = len(rec.panels) + # back out rec_diam + rec_diam = ( + tube_dict["ass_tube_per_panel"] + * tube_dict["tube_mult"] + * num_panels + * (tube_dict["od"] + tube_dict["spacing"]) + / np.pi + ) + rec_circumference = np.pi * convert_mm_to_m(rec_diam) + print(f"Rec diam = {convert_mm_to_m(rec_diam)}") + # This is number of tempsteps along tube + delta_T = 5 + flow_path_p_loss = np.zeros(len(rec.flowpaths)) + for path_key, path in rec.flowpaths.items(): + # for each flow path in rec + print(f"Flow path: {path_key}") + total_p_loss = calc_path_total_p_loss( + rec, + path, + path_key, + delta_T, + tube_dict, + manifold_dict, + fluid, + rec_circumference, + num_panels, + ) + flow_path_p_loss[int(path_key)] = total_p_loss + print(f"Total pressure loss= {(total_p_loss/1e6)} MPa") + print(f"Inlet pressure= {(outlet_p + total_p_loss/1e6)} MPa") + return flow_path_p_loss + + +def set_tube_pressure_bcs(tubes_dict, tube_pressure, pressure, cyclic_times): + """ + Set pressure bcs on internal face of tube from flow + + Args: + tubes_dict (dict): rec.panels[panel].tubes dictionary tubes from a panel + tube_pressure (list[double]): pressure values for tubes in a given panel with height + pressure (list[double]): properly sized array of pressure data to define bc + cyclic_times (list[int]): list of times across all cycles of analysis + + """ + for tube in tubes_dict.values(): + # for each tube in panel + tube_pressure *= pressure + tube_pressure_bc = receiver.PressureBC(cyclic_times, tube_pressure) + tube.set_pressure_bc(tube_pressure_bc) + + +def set_and_downsample_tube_temp_bcs( + tubes_dict, + set_init_T_to_inlet_T, + inlet_T, + cyclic_times, + num_cycles, + analysis_type, + loc, +): + """ + Handle specific details for temperature BCs in a tube. + Set the initial temperature + Find and set the appropriate temp for 2d analysis + + Args: + tubes_dict (dict): rec.panels[panel].tubes dictionary tubes from a panel + set_init_T_to_inlet_T (bool): whether to initialize tube temps as the + inlet temperature + inlet_T (double): inlet temperature for tubes + cyclic_times (list[int]): list of times across all cycles of analysis + num_cycles (int): number of cycles in analysis + analysis_type (string): "2d" or "3d" to handle turn on or off downsampling + of temperature results for simpler analysis + loc (string): "max_T" or "max_avg_T" to set where temp is sampled in 2d analysis + """ + for tube in tubes_dict.values(): + # for each tube in panel + T = tube.results["temperature"] + if set_init_T_to_inlet_T: + T[0] = inlet_T + T[-1] = inlet_T + tube.T0 = inlet_T + _, r2, r3, r4 = np.shape(T) + T = np.reshape( + np.append(np.tile(T[1:, ...].flatten(), num_cycles), T[0, ...].flatten()), + [-1, r2, r3, r4], + ) + tube.results["temperature"] = T + tube.set_times(cyclic_times) + T_3d = tube.results["temperature"] + if analysis_type == "2d": + if loc == "max_T": + T_max = 0 + T_max_h_index = -1 + for i_temp in range(1, r4 - 1): + if np.max(T_3d[:, :, :, i_temp]) > T_max: + T_max = np.max(T_3d[:, :, :, i_temp]) + T_max_h_index = i_temp + tube.make_2D(tube.h / (r4 - 1) * T_max_h_index) + elif loc == "max_avg_T": + T_3d_avg = np.average(T_3d, axis=0) + T_max = 0 + T_max_h_index = -1 + for i_temp in range(1, r4 - 1): + if np.max(T_3d_avg[:, :, i_temp]) > T_max: + T_max = np.max(T_3d_avg[:, :, i_temp]) + T_max_h_index = i_temp + tube.make_2D(tube.h / (r4 - 1) * T_max_h_index) + tube.results["temperature"] = T_3d[:, :, :, T_max_h_index] + + +def process_single_panel_analysis(rec, struct_output_dict, single_panel_analysis_id): + """ + Remove all panels that will not be analyzed + + Args: + rec (Receiver): object we are analyzing + struct_output_dict (dict): dictionary with output information + save_struct_to_vtu (bool): whether or not to write files to vtu + st_fname (string): filename for structure + tube_fname (string): filename for tube + single_panel_analysis_is (string): name of panel to analyze + """ + # file name mods + st_fname = ( + f"sing_panel_{single_panel_analysis_id}_{struct_output_dict['st_filename']}" + ) + struct_output_dict["st_filename"] = st_fname + tube_fname = ( + f"sing_panel_{single_panel_analysis_id}_{struct_output_dict['tube_filename']}" + ) + struct_output_dict["tube_filename"] = tube_fname + # update model to use just one panel + single_panel_model = rec.panels[single_panel_analysis_id] + rec.panels.clear() + rec.add_panel(single_panel_model) + + +def save_structural_results(rec_struct, st_fname, tube_fname): + """ + Save results from analysis to files + + Args: + rec_struct (Receiver): receiver object we analyzed + st_fname (string): hdf5 filename for structure + tube_fname (string): base filename for tube vtus + """ + # post process + rec_struct.save(st_fname + ".hdf5") + for panel_key, panel in rec_struct.panels.items(): + for tube_key, tube in panel.tubes.items(): + tube.write_vtk(f"{tube_fname}_{panel_key}_{tube_key}") + + +def run_struct_analysis( + rec_filename, + set_init_T_to_inlet_T, + loc, + analysis_type, + is_single_panel_analysis, + single_panel_analysis_id, + inlet_p, + outlet_p, + num_cycles, + solver, + struct_output_dict, +): + """ + Runs the structural analysis on the receiver with a given set of thermal + hydraulics results. + + Args: + rec_filename (string): filename of the receiver to analyze + set_init_T_to_inlet_T (bool): whether to initialize tube temps as the + inlet temperature + loc (string): "max_T" or "max_avg_T" used to select location to downsample + analysis domain if lower dimensional analysis is done + analysis_type (string): "2d" or "3d" used to decide whether full analysis + or downsampled analysis is used + is_single_panel_analysis (bool): whether a single panel is analyzed or the + full receiver + single_panel_analysis_id (string): name of single panel to analyze if + is_single_panel_analysis == True + inlet_p (list[double]): list of inlet pressures for each flowpath. This + list will be updated based on results from this function (Pa) + outlet_p (double): pressure at flowpath outlets (Pa) + num_cycles (double): number of times to repeat load cycles for the analysis + solver (managers.SolutionManager): system solver for the receiver + struct_output_dict (dict): dictionary with info about output for st files + save_to_vtu (bool): save files to vtu output + st_fname (string): hdf5 filename for structural output + tube_fname (string): filename for vtk output of tube analysis + + ToDos: + rename loc to better name, make it and analysis_type bools? + """ + # passing filename here because we may change receiver significantly + # i.e. removing many panels + rec_struct = receiver.Receiver.load(rec_filename + ".hdf5") + rec_struct.days *= num_cycles + inlet_T = rec_struct.flowpaths["0"]["inlet_temp"][0] + # times we actually analyze receiver for thm + analysis_times = rec_struct.panels["0"].tubes["0"].times[1:] + pressure = np.ones(len(analysis_times)) + # create cyclic time for life + cyclic_times = np.tile(analysis_times, num_cycles) + pressure = np.tile(pressure, num_cycles) + for i_cyc in range(num_cycles): + cyclic_times[ + i_cyc * len(analysis_times) : (i_cyc + 1) * len(analysis_times) + ] += (i_cyc * rec_struct.period) + # initial steps for pressure and time + cyclic_times = np.append([0], cyclic_times) + pressure = np.append(0, pressure) + + # set tube pressures and temperatures + for path_key, path in rec_struct.flowpaths.items(): + # for each flowpath in the model + tube_pressures = np.linspace( + inlet_p[int(path_key)], outlet_p, len(path["panels"]) + 1 + )[:-1] + for i_panel, panel_key in enumerate(path["panels"]): + # for each panel in flowpath + set_tube_pressure_bcs( + rec_struct.panels[panel_key].tubes, + tube_pressures[i_panel], + pressure, + cyclic_times, + ) + set_and_downsample_tube_temp_bcs( + rec_struct.panels[panel_key].tubes, + set_init_T_to_inlet_T, + inlet_T, + cyclic_times, + num_cycles, + analysis_type, + loc, + ) + # remove unnecessary panels if needed + if is_single_panel_analysis: + process_single_panel_analysis( + rec_struct, struct_output_dict, single_panel_analysis_id + ) + solver.receiver = rec_struct + # run structural problem + solver.solve_structural() + if struct_output_dict["save_to_vtu"]: + save_structural_results( + rec_struct, struct_output_dict["st_fname"], struct_output_dict["tube_fname"] + ) + + +def create_receiver(tube_dict, num_days, times, period, panel_k, num_panels, results): + """ + Create an srlife.Receiver object based on user defined inputs + + Args: + tube_dict (dict): Python dictionary that holds useful tube definitions + "od", "t", "h", "nr", "nt", "nz", "spacing", "T0", "tube_k", + "eta", "tube_mult", "ass_tube_per_panel" + num_days (int): number of days to analyze + times (list[int]): list of index of hours to analyze + period (int): number of hours to analyze in a day + panel_k (double): panel stiffness used in SystemSolver object + num_panels (int): number of panels + results (list): set of results for tube objects + """ + # array of 1's with leading and trailing zeros multiplied by input pressure + rec = receiver.Receiver(period, num_days, panel_k) + for _ in range(num_panels): + rec.add_panel(create_panel(tube_dict, times, results)) + return rec + + +def create_panel(tube_dict, times, results): + """ + Create an srlife.Panel object to add to a receiver + + Args: + tube_dict (dictionary): Python dictionary that holds useful tube definitions + times (list[int]): list of index of hours to analyze + results (list): list of results for tube objects + + Returns: + panel (srlife.Panel): panel object + """ + panel = receiver.Panel(tube_dict["tube_k"]) + for _ in range(tube_dict["ass_tube_per_panel"]): + panel.add_tube(create_tube(tube_dict, times, results)) + return panel + + +def create_tube(tube_dict, times, results): + """ + Create an srlife.Tube object to add to panel + + Args: + tube_dict (dictionary): Python dictionary that holds useful tube definitions + times (list[int]): list of index of hours to analyze + results (list): list of results for tube object + + Returns: + tube (srlife.Tube): tube object + """ + tube = receiver.Tube( + 0.5 * tube_dict["od"], + tube_dict["t"], + tube_dict["h"], + tube_dict["nr"], + tube_dict["nt"], + tube_dict["nz"], + tube_dict["T0"], + ) + tube.set_times(times) + tube.multiplier_val = tube_dict["tube_mult"] + tube.T0 = tube_dict["T0"] + # initialize pressure data + pressure_data = np.ones(len(times)) * 1.0 + pressure_data[0] = 0.0 + pressure_data[-1] = 0.0 + pressure = receiver.PressureBC(times, pressure_data) + tube.set_pressure_bc(pressure) + for res in results: + tube.add_results( + res, + np.zeros((len(times), tube_dict["nr"], tube_dict["nt"], tube_dict["nz"])), + ) + return tube + + +def sample_parameters(num_threads, verbose, rtol, atol): + """ + Create a set of srlife.SolverParams to be used by srlife solvers + + Args: + num_threads (int): number of threads for analysis + verbose (bool): whether solvers should print verbose info + rtol (double): relative tol for solver convergence + atol (double): abs tol for solver convergence + Returns: None + """ + params = solverparams.ParameterSet() + + params["nthreads"] = num_threads + params["progress_bars"] = True + # If true store results on disk (slower, but less memory) + params["page_results"] = False + + params["thermal"]["miter"] = 200 + params["thermal"]["verbose"] = verbose + params["thermal"]["steady"] = True + params["thermal"]["substep"] = 2 # 10 + + params["thermal"]["solid"]["rtol"] = rtol + params["thermal"]["solid"]["atol"] = atol + params["thermal"]["solid"]["miter"] = 200 + params["thermal"]["solid"]["verbose"] = verbose + params["thermal"]["solid"]["substep"] = 2 # 10 + + params["thermal"]["fluid"]["rtol"] = rtol + params["thermal"]["fluid"]["atol"] = atol + params["thermal"]["fluid"]["miter"] = 200 + params["thermal"]["fluid"]["verbose"] = verbose + params["thermal"]["fluid"]["substep"] = 2 # 100 + + params["structural"]["rtol"] = rtol + params["structural"]["atol"] = atol + params["structural"]["miter"] = 50 + params["structural"]["verbose"] = verbose + + params["system"]["rtol"] = rtol + params["system"]["atol"] = atol + params["system"]["miter"] = 10 + params["system"]["verbose"] = verbose + + # If true store results on disk (slower, but less memory) + params["page_results"] = False + + return params diff --git a/srlife/receiver.py b/srlife/receiver.py index fc43b01..7dee2c2 100644 --- a/srlife/receiver.py +++ b/srlife/receiver.py @@ -1509,7 +1509,7 @@ def __init__(self, radius, height, nz, fluid_T, film): self.fluid_T = fluid_T self.film = film - if fluid_T.shape != (nz,) or film.shape != (nz): + if fluid_T.shape != (nz,) or film.shape != (nz,): raise ValueError( "Film coefficient and fluid temperature data must have size (nz,)" ) diff --git a/test/Rec_Gold.hdf5 b/test/Rec_Gold.hdf5 new file mode 100644 index 0000000000000000000000000000000000000000..1fb9fc4f1af9d0ef7b9ef4e7c04e2200a4311692 GIT binary patch literal 115736 zcmeHQc_0}vCg_6*PmBk@Q(*VYviwR-(T*GJesi;- zLwTSi=Es)(nwKe#K1tvmLffc3U^GbJ#dzH+k+hZgwu7Jf6MFEFX6p zFGui^`NQ$4iwpiuE62tw8P`%gO_sMAPl_kObGP$wnl*#RlLY~aVw7gTM+d^ytj38k zxEVHnOcm?G;d)E)j99)FJb7>gPcJ8XdpA2z+`jfAB1}!;)r&Bw78uW5;Aj`~y$JXZ zzr)uhZl{tcsxQr6Z}nZ5SshUMuX_!&Px$}W_W*fa;&`%V;|F^HT`fFbUv>wMOBxpt z(2g5pVkp9E@@2GIl!O>DZz8*$lY;iJq0g&g#~oabjCTR6@zWO{KEz}c;|ulz#(qIJ z1U%jhC(r4#+}!M@d!eS{eoTyQxrdX36X^DMtXshuL4W-Wk!8CRwol~ojO=W@hPinA zVtjkH{4ggM*CDg)?CqVVJK4E;4VmRO2OQIBjvf0w4kwbo8T&mxznC9y|6Ts(EDx?1 zZWs0gd0fdljw{9P|GUSPX3LZOGIj8`6p$bH%lN$F>ram5!Ta&`$?oQrlP5j;yU+*+7AfQG4jmVex+3sxIVqP+i1*C5@#YMpE zC-fOUZ}`4rQFr75S@wl^ef>aOFI4_Z|KDglC~PS70PY8pS$DXz)`I;Y!V|AG&P;t# zL0R^~t!rG$nUBgnUdR%&hP?0aRYG6jejtwB&gzQx%Ob;5TKaW4BX-+yIa{ES8GPO6 zEI|#oVB0Un(`7!vIbN9W;KA;fQ=-_um*#0p^0Ys*6OY*M7tVvfQ9Wo9&rVcd~( z>OQE&^$MywPooKwr{TyRzpgxuf18gDi~< zFWDGz2V|z;QKKJt=eoR^veSfaSpUDc%|y2dtX>qSAxF zs_IWQxDG6A|FY}r^aZuaW%t{T88cO1P-}9@WX2iZRDD4$tRMU%{fzeVcs?0cmIoA5 zDf4a~24_Q4sQJ#?a;E%biq_aJkv}(?@{I}nyrgRih09}mG`mT*F~Dh`x+P-?aFIF^ z14@IEDaGih5Fdb>UqSn~uLW22o3Z$Z&Z(4s+Yx^4*QQVrDNg${1CuGq$>!?A{gSEo z$_tw2E7d%nGS5`)V#ic!r^4j5A6BPO{T{T4J-8y7x_xy1vync@RNlNf>o0VviQD_A zWZW@@R7!DstJ%#~rcfhVI`Jx&CR1^Jd`p&lB~!Tp!PZYY*2F!Vl(6Nhd@5zz`|iL2 zfhm-F;JGGNi;}5cuUv;zx+hb{lGGd>g_^h(Epv?|WmBmePgiz}Seim@k}$aZd0sL# zd#hou7?)(~Mqiy4Nv&(*T5hmeY}+E058OTd-+>C~%H61m}>(y2*Dr2@SdWz;x-X+1Kt zF1E;|KA5gn>J4zl>hW?T0M4&xuU%3AXZ7&Im52T{ar<-TZ*SE+lRB;^Zo6tiI>qZ* zl(A!OI(4eQ>=3UG=~UUk6GyZb*2L*eySAo}cqUckpT8!5TskGb=-R4EpLEJ0)p6Ik zcInhH4gG0W^J?OPeJrnuF8CQ$MDw~RCHOi z=Cq;d)P_#u4^MMSr;_yYO+1>XQ!jLbWEEX%;_mkJl40 zor-bob&}sCow};}ct*~QnmBwt3GL720X#klI&R$IFAhmkcA!p3>7`@=+vUyu%a?(z zFq-qS5Li=@s#dRn)mC)qC#K>+?b>`|Tw7qVck1fL(f8$F%GFn%CjZqy#^C|{zQKZRF23(IW$#P|!yXc_p?|YD`5rcr@c-YWY4pYF z>KDo9w)*--Uqj?d9CwZ0v5m|DWWT7yw%?KM9k4qb9A|lYS&_YC_46W89K-6gb;U3F zUW6gY@}PwlctGd}VI05<;RnWn#X`p?eZYGxXa7wby)OjOeMfa3uicDKZ#uuy``~0g zoipXl60-z8E$VuH=v&EP`pAHHuU3HfWk%k949??|Y4ccsyWsQt@<{;qaj|Dm0l;ar z)=ZoOa6W@3j8>Eirb{ZDS3=xncOP^-lR=_u0M1vzc;`rfyWzDmEf3(fb*k#+1aN+x zPq&i;xHa+;3~tbH7Y6t6$ZLH^fLk^#@O3YM8xvqDp9FC65*?j~0o?69V*{UpdW!G= z$kg+=)f5KT*Y2pUCcrI9n=<*8IG@h{^I7CPfIGA3wU#`eEu*7 zcSv@{!;rs-`47yP*&gX?O0q#zdH zRzw2T!$MZkK;GWJ!xPxgfIsO1=+TG*OSb#ex zm1lGl;1o<$RL!>64)?4H`pFql5LuDKsP4;xpzGPv!>vPUEUT;||T5f1^*R>6Gvbb#x= zI>j;x;NIwj>SzMo(c2vuT=wvp46Z^nz9byr%metbR{(BjPu(yxfJ@u$^H(guE$#R6 zQWt=mdq{!71-`IkaMXRivjf0gQ!LrN1K@(@*rY21T*~AvAH4x?n#t0e6+?n)^AXYv z&Ue@-24`cl<9rW*ds(d7&lTV<+?g=xX$Xif^9~*L0IpBl&c*uxPWsJzU)a7o8~QW2 zRU0~wx&q>)`RoZje*?H1M|0fc0B)vvf%J>pe7Zbdp63W~FUlSV zuK#`Q_Lm0o>7tKnE5-m^Mx^a;n*mO;iTab)0GG7v3WMwMOp(FOe=w}f1mMJ|{$&XO zC+~jKsu#dj9ZW8s3vkkhEK{Bs1k?P}c?>SJ!#w6a_im5vB8S0w4E^vVA{W@vN1eAm z1vX-(zo|HQ4xhHaLS7kI&pDpbH^BbmJrJ(d(#5SIF0sd=c!={UZ?_lXvXcs}Ax`zc zxSkL{fhRF5Kx3aqHIa?E-O9rUnHN_hu@>wP=cP>m2T!g}A6@N(}D$h}jJ8 z;)*2#Gp0#DU~_4NYscbdS?8!j z+`Y-V?;);9f^vI^8(=>&2I8U@YVjcMHH*_2f{ypB=(Q5WwU-_D9OAT2XiGz!xcZva z5NEOCR0XW3&H!}0?z#xqA!CU&#A&QB%!jyrCNl3|1KZ>5!TAuExa~m!#EsQJI6YOi zeFtAuW`lB57TIM(+@Ac_MG)8K_g)SVcWYwvvk61+w#9bc# zAs6CCmAQ|DxRe1J`yeh(0pZfSA)I}${Cg0mee78(#0}94VQ?|6WEfm`-`u?rH?1SW z{oUtxQ8+GxKm0I_f0v4%-|POK>+`$gNOZ)k-;l-sf9&(S;}CpYp%ugf^*xt6jzDZm zf{+T!{xTjQpUc^4caS*sbrbDfyq#v$@?4Jg&#`lvHQmX}ulhM%`D=KMt(yHELWrPr zj4$~cRxP`_IeCmH9ow3>z;S_zer^(O(z)K3hNVd7G)Q@ObDn3jHAT z1H2%FaUhHXVH^mr17SZP><5JXz|YqYz;pJU{yM76zNi7-+a+lyK1v7*qA$u^9T1rl zL?^B9sMWg@pKh(_D?iwoPxqbSWU}ELpWYp{?^ytN-%lQX@P>s9Pq*#uC*%O0`^B4E z>;*iB`$-lL2RyHeCpDV_cm`?cW*-1N;VqBX%9GP%0mGvcJ#X+Kz+-!D#-l{QbHDG7 zoPmJHIFP4d3V8M%&7P73cr-$@YUMFd@nU$EeAaBY5%6R#tI}HucvRnZ?xzHJ?%(e? zRvYki{AeN;3U~~rW!B1LecYbm5i7HkaR)ph=M*m30-l5mY7bs_3?Oj&ilY1x_O@qL*HZv(Z?&Y5`qDb)hc=at$=5A?1EKq&+_T$_p58= zS@Ij=`E96bw?e?PdcR4J0>ERN_H>OG;CVP~)$5gjXW;&#-1~sXE3$STG7h+IbDfXF z$HVig*a|woO3t%#oo^-QUAfM`lJl@!=VQrvS+4W5uJh64yfoMOX>y*L>wGmiZ_Rc7nw-bx zI-gC>YjcesJkL$Wk)ZK|@!#9xc+RrVd!05AjK6$?=+CZT98WElo_mZ>%T3eTw-w}D z+}9o4D#5!te}b$Bj^{I#?yvFzPwRkEx~ChsPg!ohJ`V8gQC>Id4B#0RD5}1$3Isejj!7OI&+{{P>l*_e@8?P)v4E#t(YVzL zfTyqF8NMOl@vH1pZUA_29Fsgao_!_`z426!PrHuNyg3~3sLhHGJf6j;b7NPnFK-b{ zhjy3TTIyAuk3sU_cwVVEY4Hlcvt(W`ml78~U6%1qZwBD`6dM&04tQL8X5|F|9vsIc z502-Le;!d71$gcy{ppqlc&;gLH|PU+=1Obaw*owdJN8XB13Y9LlY9+jTV|RIlaFBA zyUT%ia36j1-lv=J_lj<%VkfTyn-{R-^b26;MBFrN1?ItT^^Glt`^oyhvO_%MrDx_q zo+lwgiXqRsE0&jG`Jp#XUx7SBj~+LJJW&H`=Ml3tB0) z6$vxgv$S+vst$n!?McAkSKh{wUPT`c4ov_&)v z@<=B?-vN2VPRi|uJby`NsX(4%-D~HuGeSJ==1*S@dFEf3lLC3dPMNKTJQ3=LVj<5g z#jp;LC#rMpJd=zOPf}YYhR6T99>cRDsWZdV({dHVqw~@M@hH~LL&gEuZLaw=_;@%k z3R^+*qsTlduK7}A-W1pTDKd|WYd#g3SH(5Iip;a(nr}tsU2)C7BJ;4g#t+WRBI8KV z_`!H=x#kyPJkArtR?vJgGH;A){ur4@#xDpOJZJT=UV$ zyfm)ygY(qLI1)5|F#da6uJKIPgKIpKJY3_Mo=G0A@l5jk z)9wE*pDD3H<}=locW4b4`2>Fw$O!*8QV-z#Di1qPu)>rlZ??D14D_>TE4+*GjE%38 zjm!V-n&s_f=V4_#9lOm5Xb}DXwmU$3f*)7oBL6gR3$K4lu76a$PL>$6zLl4sJM)9? z`qs5#;+_qDn6pi;M<)1sW+-zLu7?Ko<{Vkakn3JyKh(J%*}uH5Sj~PCW!`|8r$xgaiz<+K%^;(c>Mba?7!+f zV&tis`?=;9JF@Xdi+z0#@D=&>-hXSWzVX+cT?Y`i``@aEkhfuapw9efcKttGYn0z_ z&a(&ldvG5fKWJP)K#Rvwa=k&W<0zw}#|AOH;JA*XW|>`+A6l_fyz51?y76b0GLSEQuQERH1&c`m_aBe{2As$nh|9-nT1-`!8$o+Y{kmwFiD& zTowEjaaHelxCr?U$5XW@g`XOO2kIIRCH%i`J$_#OUcU#b?>fczt<_f0^^6tJfc8_Z zXI!n@oVoaG|C3wzdbK9b6IrgYZ3`Y(T<xi$0bzgz+JL8m4B;a-TtM3<^GVgE9 z+1Au`e}(UV85Zz6fx2`m3YMee{9T<2UJxpe;~8xU|2Gy7;PEG!ZLm9QZ z-9+Elr&FPD7nRd3KozJ&#uFD^$gz`r$ZG+Jv_0&B&Uv)emI*e_T%D zI{HGHydG=g=>USQpB6u_x_q&?XPaPrUhGq{1>+!@@Q7jidWbPuHu zHn-{d_DmQ(@tptrjsTaZJ!$13P|u*i>03krPA0;S!O3cV+78<{Az62xD!}>sPf>0I zaQ)w$9W{P-1f8JNEupl-I{JfOkYYEnXnIFzYX;}8?NntKHOpxs>iqW0^6-6pf}I=dQPae#ugMy$v33^{LiRHM)k<^NyVhy>dE? zUK?p9GG$f-turyq{kp%N2e`xK zLuEd?MbMI`B<{@xIQ{g#nd<=V#a~7YZhRk?aCkhsDI50;y&Xy`Uk^>rTo_L0Y*=}D zrg9`LqLQ{#>(V;+GbQrjR!A+RioWZ45ESS~zRwyk~eq{5C1>yAkuEyS`y&~zI8wakDySR=H7i;#| zGA)|!_*9m`HI*$K-Tgtko zM$_IYYhOaYJY94@6ykcfKJ=pEdMMqpWU^w_+;DoR!NKQEdqmQuPSdaMKDUlmn0+iI zKRKHAXsW{C4$jPHaG{R}4V(jT8XGhAdIQ{!cNZot?H);+(?d^c0Nnay^piUJq*(a~t4x>kU3-*Uv9rjZIxG3lO2F*ZuEs=v=OMxA`2I0ch5H4q= z%UX!j$#k6!ai!_o+#$~0b?F}vr`Y^$dx*=5K)Bb>P(9^u&+QCx)^-M;U_FQZULpo@ zhb>}qVLja|b@sw~s#zeMMgYQzP3!Im+qbzqb_B#3>CLr*?K`+qTpi+;DTKFyIQ~l+ zrhUe82-jqU!P=VuwQ6Q7U0TMGMvRgdC`#0lO6-Rb|J#8+Mka6S*m{f z2Ed6FDrHFn+|upN3M+`8)sC+eJ1264;Yz5NT~jvw*i z!Q-uVD!l@6B96OHL)^#_2zPAN8|HX7MklOS1v)$#&q}+_061~^V-+p{_q)AY!emfh zd7f6Jo)*xr+w6W1p1-D&RSeEm)!ji2;P^d7%OI|m>6P8C02jMxg5?x|^URAIJ`mvY zI?ZNqt7H+*?M5$4h$~ov5chti&O=zwR;rShVLe+Fpn7hcjBwAwY!xBS&H7g2ba=d#hW@a9BmTOo z4RM>D@AZT@Yj=bzy78K6U$bQrfAog!yHa_=7T}(ETJM?(aBW*Go-qmFcDkfZ9RP5J z<+GUf_4Y-$CGpMXr~q7x6N8&Voc>g!*UkW^bIBpq0^qt>oJslv;O>VZ+>{TX5*G@13KUBRCyiNIE#|{!6z(QMA0zgzU=*!q?^7#}X{;#2)2@lZ-A`lJ(rvF*I*S4x zAF21!ih!qi<$Mo6z|(!a)^RwhUE!X*4a^9Bf{4F_;%Qb%Rye=6>g2oTVW6O1Z665haGq!@x zH6|uz z=R{=8-gbaTJuqDV6{yd&PLlKV0grje=le#02gfnVgX1}6H?8gp+U2r(x}`a2SHRIo zjUj+%)zz@OL7-i;EmV&Lf_C9JCV6l?XTLQ2({FAVUHOTBMz?bWJ#%Bq9fLEGv{$R) z!~I>OXiwSm>h5D=zT{t!JUE`m_{-}W`GwKrZ&x(E)-i&faG=UvBqx&2yE6D-UzaF) z@|p7WtwzUu$;Tjha6E6{n>Lu|8%8_pmpe@B5JC5!7-XAxI+9iiR=#x6If^!v{L^gP z$e1tr86*#mXW6{>CnorW(SPW;NbP7JLHGG%MX$ho62M+SU=whr8$JDNO?N;q4?uWN0uSeju&b zQXT9k+j0CU;(7ad)5QLOr{CDpOJZJT=UV$yfm)ygY(qLI1)5| zF#da6uJKIPgKIpKJY3_Mo=G0A@l5jk)9sqSOH*bj|9u_p zxO1rPi9{28iR~Y|iPw+9zh8|z@$KU)v0qo6`gXa-_U(TcS4QWNuNzie)tl${7pgUm zyT%qTjlI2Oo*({QvJv|`e(ay`y?zwQix*99=<6s;5bGRNuZO{X9c8#i;y+)v1Z8M8 z^mUZ+dG81Rzn|5d1e7#+T)T8&? z`WsmnQFdKUpOf8NL@B6W&&!=tL`g`Flka-2kV2)a|4*nA<>%}yq=GDaM4hxMqPD$r z>>Zp|L=`WI-Z={3-g`z*+jF)iuITip%vn1MsS$?nPyJz8MD6W-bkw_)B5HS;o^F;` z5v8@^YkS{rHFEVAUn)5sfbdy$z84LQABAr(>f!YT@&}L z<=azI+X^Y2;Nf>!pf9u$n*ey~?5(tv`BoDdpOA zrg3EVQfjco1&^h!HE~;48qF#6FQEq9?exJ;wUoNYJhdjyzij`ZEi+4~3Jq!5*tVtA#P|97 z$Gl3ZfHz$~&XO#p#_J>ukse&9?!__{PL*rU( zxg|5M?Q}rnn&?C{u2~I3)YiT+qX5x!>{kN z;JYXIhb?!Vj^v*ruIhUqo^>Dj1IH8E68>*29;jlbV0A>_~bJS6epYyXp5_WHs~bzgz?6pWUlxNv0Gr>?$V z2fq)x@#t#@czN5}K`CPV z{?6#@uLa)1*F_Yor;y1WPYSO~#cRRnjCguE+1tATjav2y_36cm!E$UcHwPcY8LEqX2$!6D6y zV&fL`l5s2`eWNK}m;O-VbL4`*t3QF4S>CF-9|PxWq$GBWBjXsqdAWO657O}W59X4vuhD)S-E5Asr{>$7YO< z(-<1*ICY5W&geL8U3-whEn;ahe9VeY+q~oN29_cuV z>_9qByG9}%r&pp#$EjQo={SvMb)3eHMmkQ{6OfKm>+wj(DI_VC(QyjDCWUmIh9MoN z!K{wcg-J-qNhJj7I4z1mI!` zCuMD<Nu?|LOM=+vXG9`fmD6gNM#sr?iyx!owC%y~3~nZ?Hqvo2u|PUbZ>*7yQ~SyTjE>XYMFNvTM zKsrwI(vXhR_9Ud^bbLC}ajN=I#ppQc1t44{tK+o$2GVg_Z;f=E+G-*lCkYdz<78Tb zbez&hAe;`XigcV5-B~>*R>$dXG176;$wxX)e08MbwD2_2amp!1I!+-S5iWq$ zaSC-nI!*^Fk&csV2GVgVu|hgdB~y@&(}e*Dm%-{dEggw;oIX89I!>pgkdD*psYu63 zeHPMj+7^m*jw)LqJtkQdq~mn%H>Bfq(Fo}{eVBuEoF<({I!@LyNRMePtK%ep4e2;3 zFGM;{aYK-fQ^;MUgwQZ&+Wnw{`%9j9F4y+M7?<-X0(n*{vc;BVaZ3XA+g9w48~ zmDsk)isLs-Xh-AGaY|(0?ESYoPL72Cf1!~3SDRAT`%3bCZ@v9u9~A#e|7vZ8SIWIL+vQbe!DJUSxEfDho=Hj?>IDwetMQ>NpK~f^?i@b&!tJvg;*`j+4Ws zJ&cZ%*~31Jj?>dQ*|qY_Vs)G*zW;;Kang@OI!oaT!n9j7WAq~kO*7U?*JKAywqIK5kE$>=y8IM%<``V3@soQzH(9j6iYNXJQI zG176mlNiV7IB69YFgi|BchA?Vo3xwNak{C7bev>Ns(okHeP%o>#?I z(D~I8VxE=jd@DKc%60ygoQLH)A4|^5a-E+g=V`gl*OK$LT<34ed0ejXgXeY0I1)5| zFdkd3^OG2l=b5n;biO%(n0Mwn|4hz9bDfVS=cPH%PZRUhoad{Fd26oo*W^4l*ZFL6 zUYl$D;CXH`js%S#jQ`#i$Mc3Dq~jEnfpnbuIv^bz5rW3nDNo*%xt z%jh^gcSbr+iLH^2)1rF{jE>W?W=O{g$1%x+<2mI7(s6QdKsrv|G}3XpzoiAE<5Y0> z6QkpV zJkoJ$9zBfFagwz<$LKiW{0ou?$8+#(q~lbz3+Xtm@<2LH&XbXjlc)mHal&y-^5A&B z8;^9Hq+1~!r-dz%j?-Z|q~r8_(k@2F3CA(XgX4M9Nu=X+T^#8+dA?ZA=s2l1Lpn|) zon;suCo+ynzJ{``#OgSmo`-atLR65B)3TRH$0_>?(s5Gwz+-fr`faaW&nJ`Baf29)y{K=)p4@8j&z(>Ya$(|i^WLCXzN)rn5Ru$8?a6)4(L8 z$ekAf)3YlZ8hm`57lo~$`B7w^6xVzyGH;4&{uG%<#WkOb%&X#> zUq$9wam}|P^RCu$%)cV@u(-w#&dVa>NYMDfcx<`m7hyck6T?=}d@(X_jBEZFnMcMo zpNz~aGal-Aa=8XVG!WCDu9My6!Y{y|EXaByh$*U$+EhX#OlZPVoG5(MO6Ifb@|B zwLP9AeI!9Gj#HsNl29K>sE>rN1|CnSk5ugo;itMi@HQmu(MS~sYKf0uQ7N!f=XAIm z)Z2klS?-yutm;4&$1FaS4Q#7D`djSP9H_Xx$&1C*9lq!#;ins{tHf0w3FFcXC#-i_ z0Q3VV_%)5c97v}+Y?yhkM-Z*$cl3(F+#ouYU#Mqx5o~V|+ZNoF|3&}rgiXufF+g{4 zT*#9`Ctx2O$sVvASkGr$L*syL+G$bgL14p-ry5)b7J7us__nj`e%mo)rs@l759ON7 zIK!K&FQ~<5t}Z`wX*c8~`Iv3?eh04QRP#9v-FbdL9VaDL;RxT+;+BTs4-LxV=YQxM z_?Jat529`IKYm_w9Q5sSB<@$=t{GRlOp|}Sj+3koIx`Jh$0=|qswIx6Up+pA?fcpT zbw2;I>zv}!|CWvuxxOIR=js2~agqX+=K4JS>+3k_vg0D^7xDRR=sHfh!@fPv;r`3l z@Z068+uJxk)f^AMzR$w?P1thRaWd6Io^RN3Rp0yYz>&i7uyI@&wYRda@lc6a$NPJ~ z53lb!{pUJP?mzYOQ0}O2_v>H!f8spVewG{0UcuvvtBzA(9q|>|pq!2MYRGl=tM3=N z>NrU;?>x}x$^$BE=t;hcU7!E&>QX3}3&)=Z6)-=zGI;z+W{>O6+RShi7j@}UByU6C z%d*{MUF!fT??45NcCiInejKmLNoc zJ%vE_c>(vBY{2nC@P`IvNj(KwXA}WA4s_XV3m3Gbro z4wq93MW1N)y0>P_)~eAeT?E*HpdCr(8BaePaRL)$Ej zv#G0JB9B|&bDHL)uj|LTLPFqA=>a@W;&BSceGVJ9xKWKqPr-LGx+4BtJq0C;@6%7< z^Vg`XaQzyZIp0D*2>k#9;RnXS_=S#7`hdKR>HeEGdS3{ldyMKlUb`8eUU7b<_rb}0 zdef9QOUx4Zbe_xkp>HLF>EHqHUabK6Ax7SQ49@YBY4ccsyX^D(@<{-9Z?R`k0l-PO z)=ZoOaOQ(1j8>EiruS4fuY|Z&};HE1W?;Htm$Gldi)r>b5~0B6$q zbUQhKbC8!{aB9O{7~It(uk{@P&ShHQ>s|n-5MU{v1aLDYIyw&nxU4;61D}F=Ud#N* z)bp^_6b9GH?x?ONz`3POnfywePapmBS>!x`+q~$tmOQ}OSk{=zd1u6!=SC8W9J_ygSV?jDE60^D(_ zJfoWcS3aiXt~bDS6CWTW2XNg!=Pa^QYo;?%_NMJ{`FI*FKplC4gDFM z!-meIu7Egs>OP_8Zvc1vXpVavz!{krNWZwvr}xCm^Be*0RM{g2=f78p!O3ODkDUu} zO82$fUm6I?eNR!49@Yvurd>X`!u+JSpvYx zxZkwu1#lM*CKt~IxCaL(>vQzfIoTMGgZS-sEXSF0fRN^VX-p zZjACb74IEL&(5ikR|b|J=qY^z?9b|faAGSKw}!Yl)kX0Tx9DxVy%2ZoV4*d{8Jro{ z6XIrtA>36K=kBsp4&pw!W7Z*hAh4<265** z+=zy_bqA*>K%DwM7jua7*Y43C;=ER&7kxIs8{%51KU@NFZP!F-LR_qwdPOf_oBJYMG!Nm_ z_&v`-oJ32T=@6&zWT!jC?Mm3&1>%}n85BUA#4Lo9Yl?8oo$s7w>M5{4!YK_B)nsr(Mtp?0&nmhdAKVUDIA)LH1!cCo?qYiO5 ztaaZ*+_NO*_7JD)F){|?Hm}m+L7dn)gwq&;aIYV~R)V+=ZO1)_IJNWI(h$e{W6f%a zbBs7u0qdzX3gHxW5w6>TCDIT#Zlz&9#O<+=dH)(%gKG!pL)^}^2L%vkt&ebnR1vQC zkr$QOpq&2-yKIQt`KWaf#Pu1{%K_qU+c!T8adFxR*YY=n<0V)YK-^MM+cb!~W%?l( z;%wfzkAt{BM`-MWxU}8~ccdG_nf1=U2XV&P&r%`I$S{P#rFM~FaOwfMdm+wV9pUPF zE>|Mn8`Sq)p5rKdE@$px8|gptxm=0$gEjl!3A?57=qXgOZ}$FMJ%vod|M)xz?YBKp zSHDQU@2$6AT!Z>WrEecj`0X$10kU7zWZUV>_71q&KT1!*F#^S-P)`AP{s+EiV9vAf zc(l+D!nhM&cOJs)0Cohz{!kbP!Z;9K2S0m10Au(&{VAx-zNi7-+a+ozK1v7*qO)YK z4v5SNq9fLK)au=dPnUM_l^^WPr$0}3GTCsBPy0pfdlmrR_mhVoykQ~3(|&vV2|0kL zpm=kOy?`gKpJefHz>`-zso503qoARieE{%GZh5>`o}4BN7#^+Yd4mrDp3&E4JW2#S zXZr5Q83=gV2J$pa0Z-7;>?uirM>I66R-XPUUJQ@dXU%pS0Z-hrD!r9}=Ur9jeoBC6 z*Zq!TwE@rj4<=%vfJb&(X01F^kJ~dm70>Ks+yT$Ra|#!10Z-5cwFj>|^698{vXWAO zCu&{$O%8xZr+infJV*Q|Fgyo>P6Yf8cr?pg&-MX4XEnEop8z~%XA9=tbmG%T3gc&J z0UonE8*1f==MQIiG_Gvcc^?==zf`|4^i6gU9afQ*5Da(>SIPTt1w6mSE?D*UET8s% zzq(eQ<-Z}GmP1v$6#|~6`%QWj0G{z_PuF+>p2Nddy5=O-~9&og5y=zMdHd1sFE&*VHb*ZF92UYhItG&xVrb-tRMx8^#3 zP0nL;ozEucwYkO*p64dxNYMDf`0s6TJZIYHy-pj*r(5z3qCdOx={Bjw(sPgT>5{2h z`?i97i(R^dTP1i`=TDIJ!1)x%D&1e@1D^H)rF2g>KHYk``T97(vr2j0s55{^HOMXI zBH+PsO!DA(PU4r9W&@teo3kIu0iLlw_Y7A8o~2SE9{T{#>n2Sz_5&Uq$0QGq=WSb; zPYDM+3nvv!EClWP5dAJ^D&SGu_O??1;1N}QT@?s;a2%67IG*h@ck3Gi9=qpCBC&v{ z>fyN63V=t-@C@G&@Jy}jQ*Howa2%67IG)E&9D3ua9?0(*rFnBW;AuN6KJa)J$ghZ9 zwZ6PXFm2jha%-tqbv_2kgX1|+anj-yfX8)SFP9P*KAoQNPHzU_DUOYb2nRd{J+tzH z01u91k_X3g$)86QMggAdNq@Sf0iImt?FM}SkBPL#eJjAzamT*NW`Ku`W0J3-Y|Bh@ zVe%0y)Vmyr2lvr;iay+{%Eh@HF+?5(gJr(XaYxb3E4D=-iKd`N6D+)vj3r5)lK zDLXR<@;r(hQVez1LR56sGY~k81Zawuf*_#$?JOi9pF+8-0Bg3<$Z|yu}9B|#{noom| zhx4Ma6*NDJ%#-4pFGc1}am}A1^QgGyQ;~U9T=T2QJS(pGR%G54*ZeCo4~uL3;Jhp{ zjs%S#jK`L1ei6pwJTYtq%@-r{#<=E>k$GfX^U27(GOqb$WS$wl;ADQcoz360sGyeIyB`8C)q3bEs?4MD@eZO^-%h`SyM>3AI5d5J*Sv>B`vbB*I?4($=*;}L>nLNF`(#*I9#BlB%)5CQoDEH( z<~wW4nevkI~69c{jfTP>i3{U?7^S%0BRP2Ao`CF71Mq*998 zTg`5^GKCt^(ur5GG?|L)<6E-aE1Aj-2)2IGu_o@>q=YS3 zH~Q+dNNQaZ*K&i+V%rv}l<4WPmI@0~D4mEco3GAErlhMjo%M1^rv9+iF&{5e6W4T; z>UmqqOzK4E^%-dv>C~ayT>`GmPp4julE@A3luk`LDi!FxD5J*tOY4!Db+JVz^}%$t zQg47WR*#n(0dRgjd+m|}IID*rt~~UwiQAtue|xLunbdJTaobfB(kWigqKqAL)2UPa zWruimNT|?pLYm-dMqhs%L9Zb`y@O`KDs(GbT`H{~=Tg#_Y2R0exIs4SaW$rDJd;TGV z(z~o@Fmz-(RcSwZipi{WYPz(?0^^qH)Oc$CnDHJpaaHEaJBGi`prXs7HKz?tr#5sN ze|VZxI+diCZ{pEBoqC}gB&+CB6L+_#mrSeY8Po(f-aXsF>C{3yneG<0=~Rqsuao>H z>C{!#$1`$f)WqTINoapA50tH39FnB$KxN44rDOs7`osOpmw`<(o%6B~*e5NjTD<}` zt+zuzF`yG7k#u5QTVTI#WjV23Q)dE}H@&5=Pg<=1L0V>PVQ&KD{_kM;ZUF0!kbHQX z!PkU)-@hF66;QX@Au&KU8nzCM}cPheRTa-xBusO zAF8nD+x_~N{+~Ec)$MGYp9GI9t~v@KZ2wQjS73v3Jl0XL9)sd(5xXvR_5C7O9R;}J z3Xg|dkB^hcu1lYjf{s_;_`}x z5YX1OUXWudVh-H<_u|(gLm0mTy8f3cPR1{O8c1a02RGK;#?8*flgGpNQLZ+go>ulQ zvwYYBPHrxCURH1wbhHDma0Krm4)lv} z)_1>vy@D-%UTiNS!YBtpzoPHZnGoYiFxyRY5<*gJ8V-!AucEgIfdeR1{ApKlWDxc>}kc+LPd z65JnL@5{65@JqwG&VHOfZu}(vsQr0@9G~iYo`^e&`VHKVG{}oVSgC)_19jbx=@R!h z-}`>7zUvu}X4rDSk0@uaXEMSYR{uucM@XE&EwBC5x@> uHXn&=XxsJ0)jz+F;Kd;i;WL0e(Ae)I6zcFx!@AA{j}NZ*RhmDAKmG^db*@AJ literal 0 HcmV?d00001 diff --git a/test/data_6_d20_hr6.csv b/test/data_6_d20_hr6.csv new file mode 100644 index 0000000..314e01e --- /dev/null +++ b/test/data_6_d20_hr6.csv @@ -0,0 +1,32 @@ +Flux intensity plot (units kW/m2) +Max. flux,845.368290,kW/m2 +Min. flux,5.230420,kW/m2 +Ave. flux,317.202681,kW/m2 + +,Receiver circumferential position (0=North; +East) [deg] --> +Receiver vertical position [m],180.0,165.6,151.2,136.8,122.4,108.0,93.6,79.2,64.8,50.4,36.0,21.6,7.2,-7.2,-21.6,-36.0,-50.4,-64.8,-79.2,-93.6,-108.0,-122.4,-136.8,-151.2,-165.6, +11.5,5.865428,11.283991,13.880661,16.842514,22.815152,32.645430,43.186667,58.565392,66.347012,72.230932,85.402586,97.840197,109.429092,121.305551,141.017264,153.904936,151.934404,162.649725,138.853229,117.681587,74.587489,43.714894,27.364659,20.800167,6.491756, +11.0,17.658028,29.897828,34.924294,38.709276,45.014505,56.841670,73.787590,98.083153,110.704619,123.382102,152.472569,186.575658,219.479074,254.866113,308.945140,333.658735,333.745597,368.698913,315.307864,273.352880,179.640739,110.583814,70.516157,54.587778,27.069226, +10.5,36.674288,53.725361,61.776878,66.756774,72.478764,85.964130,111.261308,144.960179,163.938502,186.376298,235.371592,297.652387,357.103579,420.466061,513.749207,548.871754,543.427509,563.915575,478.943700,407.648136,272.732278,164.844613,104.091638,79.118222,59.918282, +10.1,58.385827,75.340911,86.230176,94.034261,99.839781,115.841957,150.200820,191.523915,218.579634,253.160567,321.114763,408.338732,488.118858,568.696977,680.669811,724.857991,707.606809,678.833907,585.772046,491.729620,344.654885,204.458647,127.979380,92.955093,89.964197, +9.6,72.547265,89.155024,104.661356,116.229239,123.349889,143.254884,184.851984,230.298822,267.376752,315.788281,398.582772,501.378057,590.364240,668.949451,769.061696,824.309470,794.997917,698.053072,611.227515,502.933148,360.085565,211.882003,130.134056,94.298304,104.189166, +9.1,73.377691,92.286521,116.005702,132.032528,142.310761,166.880555,211.843447,257.966925,306.702048,369.585016,462.757405,570.680349,657.203287,715.064570,778.664687,840.188562,803.299265,678.705748,599.267354,483.599009,346.786118,204.969999,125.356205,92.589900,98.652784, +8.6,68.126975,88.699430,118.705590,140.683250,156.752226,185.942657,231.035091,276.232795,336.840817,412.816882,511.412373,614.727913,685.250575,718.887212,751.434871,810.145711,772.694578,659.723400,585.398800,463.408221,333.813909,199.803179,124.830549,91.094809,91.101483, +8.1,63.940598,85.630442,116.109772,143.114612,166.228001,199.362980,244.063480,289.067585,359.612452,444.521903,542.487606,635.177725,682.745331,704.719091,734.388058,785.658757,755.393345,657.349880,579.319679,449.218987,329.253533,188.562395,114.133972,84.836870,89.174000, +7.7,65.938220,87.002191,114.267323,142.468584,170.958670,206.581963,252.037516,299.190928,375.814454,463.690297,555.999563,637.906126,671.859390,694.371924,738.305565,794.258378,769.387549,671.946760,588.653729,467.603572,344.652994,194.437526,116.424957,90.995471,93.698518, +7.2,72.749707,90.784625,114.938310,142.247289,172.041149,208.226388,255.404607,307.232328,385.451697,470.630232,555.793486,631.760997,669.048445,694.993634,748.389640,819.531479,788.565498,685.087510,601.891633,496.902157,360.697599,208.416296,128.624520,104.597274,98.665569, +6.7,75.284866,92.306701,116.572879,142.513004,170.562028,206.046851,255.191811,312.701706,389.458540,468.575243,547.455033,622.168154,672.309036,702.227319,750.151506,827.207277,787.296141,686.564217,600.968183,498.867845,354.642957,214.084027,132.532529,105.996993,99.333239, +6.2,70.965130,90.513041,117.968740,141.772216,167.594225,202.900224,254.045631,316.943235,391.111456,463.535258,537.986543,612.409394,673.077829,709.334630,740.888804,807.798666,764.351904,675.427188,582.730503,469.032813,332.308408,202.373653,120.507249,91.521556,94.799816, +5.7,65.369504,87.770128,118.570521,140.517851,165.446535,201.413718,254.932867,321.221018,393.579647,461.522886,533.271927,608.270106,672.165140,712.171492,736.292361,788.051437,755.380651,664.216106,562.734177,445.375915,320.775324,191.000050,112.753438,80.974581,90.969674, +5.3,65.242545,87.630401,117.939154,140.059561,165.567066,202.408342,258.663648,324.952294,397.262377,465.003605,535.895485,612.184177,673.197282,709.606967,745.951442,798.345822,773.845989,670.894856,571.483626,466.254142,338.391464,198.446811,125.640660,89.312089,95.177955, +4.8,70.196015,90.199699,116.108580,140.232651,167.285845,205.241671,262.334833,324.229833,398.012697,471.886049,544.303563,620.755932,671.838768,701.648864,754.459887,818.378749,785.713573,689.098141,592.352752,496.439935,352.445625,210.427393,141.827345,99.332781,104.781586, +4.3,73.220712,91.635071,114.157709,140.770533,169.512385,207.393869,260.669920,313.829163,389.590378,474.782080,554.314978,629.917527,667.109611,692.568176,748.590553,814.474930,771.477045,691.755474,595.052072,497.307045,346.584656,207.445635,140.437135,96.522706,107.539481, +3.8,70.409957,90.916623,114.300921,142.161480,170.575333,205.693454,249.149957,291.363464,367.493169,464.898938,558.319548,638.446104,669.411450,691.472540,736.063023,792.200397,744.409295,675.053970,581.309402,471.844117,329.698648,201.432119,121.422787,84.903556,98.904763, +3.3,66.780681,89.950648,115.486836,142.499273,167.874759,197.647559,227.030377,258.131103,331.951006,436.722451,545.794823,639.946223,682.582719,702.179375,734.972898,783.185103,743.191115,655.591600,564.161077,455.173329,322.085878,201.880120,117.572959,83.055553,93.472703, +2.9,69.395461,90.794585,116.332867,138.255296,159.014698,182.116590,197.061057,218.908629,286.921045,390.521892,509.117976,620.450519,688.392962,715.721760,754.817587,808.516245,776.752821,660.333729,571.716366,467.782506,332.096122,208.290729,130.425582,90.604223,102.996816, +2.4,76.096444,93.912987,115.017278,127.997929,143.141059,159.085041,162.804739,178.222814,237.049681,330.724159,447.850780,569.836090,660.615014,710.692050,779.728021,845.368290,807.663824,685.291615,590.610631,485.669678,335.146460,201.751498,131.093149,90.737506,114.298249, +1.9,76.145988,92.622307,105.890143,111.298611,120.551093,129.631886,127.064487,138.460137,185.790386,263.167141,367.977154,488.420249,590.661192,662.443853,765.339580,833.682320,789.255551,705.651209,606.208101,503.061993,344.327190,207.680376,129.159426,88.261077,107.457543, +1.4,63.130895,80.433153,85.132597,87.938759,92.506930,96.425126,92.207886,101.001451,136.037633,194.051011,278.519378,384.219591,485.234091,559.577868,674.698926,729.640393,693.747841,679.421832,588.045764,492.662784,339.485696,217.604536,119.520737,85.949469,77.003034, +0.9,41.958724,59.259016,56.733468,59.984219,62.544511,63.997906,60.940287,67.682370,91.338342,130.400483,190.666579,270.068952,353.606763,411.888101,509.713656,544.559336,528.545285,554.225394,484.146205,404.930489,279.459182,187.582891,94.763647,74.438300,43.539273, +0.5,20.874981,34.346603,29.849904,33.548026,36.352097,37.111611,35.968006,40.788031,55.310048,78.728883,116.325708,165.596506,218.557759,252.014828,310.135687,325.431913,323.981447,358.292285,317.543371,267.717581,185.578683,129.770336,62.883900,51.794553,19.053695, +-0.0,6.838274,13.739533,11.878137,14.792917,17.984004,18.461297,18.648988,21.709624,29.880444,42.480334,62.727140,88.124764,112.684268,125.552273,146.843451,150.444623,149.552792,156.682002,138.905323,112.543353,75.687404,50.695634,23.663315,19.962505,5.230420, diff --git a/test/data_6_d20_hr7.csv b/test/data_6_d20_hr7.csv new file mode 100644 index 0000000..1beb8a5 --- /dev/null +++ b/test/data_6_d20_hr7.csv @@ -0,0 +1,32 @@ +Flux intensity plot (units kW/m2) +Max. flux,983.427384,kW/m2 +Min. flux,8.922926,kW/m2 +Ave. flux,396.461019,kW/m2 + +,Receiver circumferential position (0=North; +East) [deg] --> +Receiver vertical position [m],180.0,165.6,151.2,136.8,122.4,108.0,93.6,79.2,64.8,50.4,36.0,21.6,7.2,-7.2,-21.6,-36.0,-50.4,-64.8,-79.2,-93.6,-108.0,-122.4,-136.8,-151.2,-165.6, +11.5,8.922926,15.235459,15.685769,21.794866,31.948658,45.217232,54.469739,77.368883,83.637944,99.849242,119.121781,128.287710,151.895427,170.283334,164.679445,207.063765,211.601128,202.007188,170.970138,148.781458,95.484908,56.229191,34.401579,15.069059,11.714641, +11.0,29.528849,41.431577,39.220303,46.098350,61.591009,82.280134,98.675091,135.926759,149.747750,178.183071,220.106817,249.931875,311.494423,367.480888,359.621292,466.268497,489.703524,466.166995,392.716792,350.138107,232.988908,147.912204,96.486663,47.187907,39.618773, +10.5,60.938384,76.785420,71.634903,77.705316,99.077692,126.279086,152.031612,204.699260,230.785744,274.426236,343.744364,404.745647,505.332961,604.958286,603.398440,732.057587,746.701333,698.254957,587.610869,510.333073,340.939173,219.483638,155.357787,97.907739,83.279784, +10.1,87.588346,109.179247,103.311393,110.488520,137.335142,168.011649,204.633129,270.291947,311.761649,372.172976,464.858497,558.507991,677.500809,796.079043,820.935327,902.933964,898.445635,835.036405,712.882125,607.175918,410.498907,262.293785,183.012993,154.015067,122.510062, +9.6,96.287036,127.714499,127.011989,139.574885,169.231515,200.762637,248.367286,321.850411,379.415697,456.763456,563.392400,679.621901,791.182740,895.914871,958.689898,956.132834,916.681507,847.301445,731.708069,603.305992,409.601966,254.966704,172.228053,175.158471,131.236359, +9.1,90.143800,129.039677,140.242220,160.077581,190.294213,224.047020,281.724333,357.231639,429.676482,521.138343,629.886225,750.983334,835.693442,904.219990,983.427384,939.968280,883.380317,812.272659,711.540924,575.027630,396.005912,242.918301,157.120404,163.854282,117.177716, +8.6,83.518119,122.080194,142.697478,167.944854,200.294942,241.236328,306.475531,379.959550,464.564908,564.062626,662.694268,771.383697,832.136784,870.824521,936.427630,903.316042,848.285491,778.114304,694.105988,563.496718,397.992773,255.916526,177.988589,149.415893,111.110618, +8.1,84.844018,118.463194,138.658459,165.599334,202.250879,253.897090,322.203681,392.994969,485.475046,586.946228,669.439170,758.922648,816.250320,851.797831,891.460919,882.148651,841.094116,768.139258,676.133189,549.975923,387.048359,255.158521,190.101168,160.930912,126.711105, +7.7,88.606675,121.368241,134.598874,160.105816,201.085459,261.182535,326.877128,397.053303,491.820843,593.060860,664.836097,741.967334,808.589537,860.190711,899.164891,896.083232,868.887095,802.126067,693.041056,562.966183,398.573746,251.822365,183.396783,176.402568,137.802650, +7.2,88.710764,125.769425,133.170083,157.922124,200.522821,261.621018,321.133330,393.126654,485.161978,587.240577,661.752582,737.851675,811.417253,876.622836,935.762797,923.708858,892.948306,836.065051,723.803863,591.259520,419.344867,242.410757,166.677255,168.130270,128.310612, +6.7,86.065741,124.753428,134.301412,160.309443,200.867271,255.742796,310.372835,384.517840,471.875633,574.340288,662.382007,745.739866,820.111667,881.109089,947.579867,929.113098,888.759499,832.776927,739.211230,587.313809,433.472993,244.231068,177.563916,154.157871,116.069420, +6.2,86.851176,119.416419,137.085627,164.416711,200.750613,247.559469,301.901607,376.635029,460.551694,559.939780,662.615454,756.177377,828.487583,870.598640,922.294510,905.099617,854.867439,796.788237,718.561180,558.587341,412.199248,263.844891,199.314350,161.916323,122.734418, +5.7,89.944198,117.985859,139.770124,166.781516,199.735027,242.690870,301.664359,374.567335,457.908739,550.657663,660.314397,762.249796,832.090315,862.507878,898.352879,878.265216,838.176828,783.117760,690.109151,526.153235,387.453526,260.598841,194.389202,163.884132,130.553481, +5.3,89.997151,122.085845,140.100728,166.342130,198.567216,245.294549,311.857081,381.799322,466.912137,551.666036,656.961141,760.431916,830.252376,873.655748,915.045446,887.889543,865.355925,808.681142,694.228783,549.774008,393.866204,256.784323,169.279408,147.322551,121.243026, +4.8,89.016938,125.803005,137.797867,163.431607,198.236588,253.544994,328.770256,396.809319,484.309183,563.091723,654.056414,749.780361,821.193348,888.449898,944.178123,920.186307,893.232929,835.925933,720.447840,586.505016,411.281675,240.192916,164.269971,136.802472,106.931973, +4.3,89.585605,125.020490,135.508295,161.239002,199.411261,261.879500,343.723023,410.684704,501.928364,580.298203,655.398990,738.531490,807.999827,885.633338,941.081684,926.871260,891.109574,832.541985,723.872008,593.237626,413.601192,238.403508,190.497411,159.829692,113.365312, +3.8,89.184180,119.888233,136.709522,163.201599,201.690289,264.297254,345.662312,411.045793,507.533261,594.457950,664.975344,741.203071,801.621025,870.108223,910.592801,904.258621,858.150510,807.080057,699.575834,579.324819,400.026988,257.210962,196.531796,172.245066,126.789928, +3.3,88.066721,115.749899,141.537761,168.288468,202.118211,257.171147,328.525994,390.148364,490.999016,592.495575,676.400264,761.210249,810.063869,859.240873,897.831413,875.389555,835.998638,792.386018,673.780323,561.338860,393.633684,259.267434,179.446006,153.165240,123.084747, +2.9,89.637651,118.346971,146.327505,170.610268,196.608815,239.559905,294.101691,350.298785,449.723285,562.537575,669.778693,779.674951,830.921057,874.765262,930.612874,880.033816,849.270790,812.020584,679.540443,570.250074,411.013310,258.931136,157.142770,131.240173,109.023860, +2.4,95.625073,126.155695,143.126580,162.824216,181.968365,212.401973,249.688600,300.479925,389.101990,501.243356,625.314710,762.514919,840.506703,908.864416,970.330468,914.828105,882.660814,839.002695,703.170130,575.029215,411.689479,245.738549,174.539984,150.499486,109.352259, +1.9,98.987652,125.427274,123.163849,139.814752,157.063236,177.196743,201.199731,246.809825,316.670459,414.383581,537.723478,685.399883,796.632680,900.007702,949.631834,937.812431,917.069825,860.379968,731.356639,597.231557,420.878011,260.531246,197.470836,170.412804,119.088469, +1.4,88.497300,105.834032,88.435603,104.431550,122.911758,135.594408,151.351234,190.013396,238.658210,313.131211,418.000464,550.699542,676.507731,793.646209,826.210678,894.664366,902.993490,832.705332,720.638009,599.666548,414.077457,276.754368,179.212555,153.595975,114.637919, +0.9,63.003641,73.883265,51.880555,66.051805,84.542109,92.352880,103.037852,131.823296,162.457595,211.949464,287.643763,384.811562,496.902131,595.715128,619.813791,727.392819,755.892001,682.046833,597.090584,501.748379,337.279137,237.624803,118.780378,91.264270,84.510060, +0.5,32.955955,40.059096,24.578039,34.593556,49.616879,54.264599,61.490290,79.176869,97.472717,126.205968,172.268150,228.552909,303.243562,360.032313,377.496057,466.057780,497.897431,449.949907,397.557133,339.144048,226.251746,160.985653,61.623834,35.757786,44.002549, +-0.0,11.499806,14.719394,9.367245,14.635576,24.200211,26.711428,31.423467,39.928058,50.755003,65.652015,89.589039,115.675556,150.847196,170.674352,179.282804,209.431456,220.527056,196.411236,171.605219,144.060590,90.502142,60.620801,20.447250,9.534841,13.907637, diff --git a/test/test_interface.py b/test/test_interface.py new file mode 100644 index 0000000..a545106 --- /dev/null +++ b/test/test_interface.py @@ -0,0 +1,551 @@ +import unittest +import tempfile + +import numpy as np +from scipy.interpolate import RegularGridInterpolator +import multiprocess +import subprocess +import h5py + +from srlife import ( + receiver, + thermal, + structural, + system, + damage, + library, + managers, + interface, +) + +header = b"\n\n\n\n\n\n" + +fake_data = b"""-,180,135,90,45,0,-45,-90,-135, +0,1,1,1,1,1,1,1,1 +1,1,1,1,1,1,1,1,1 +2,1,1,1,1,1,1,1,1 +3,1,1,1,1,1,1,1,1 +""" + +nr_unit = 2 +nt_unit = 4 +nz_unit = 4 +tube_dict_unit = { + "od": 21.3, + "t": 1.25, + "h": 3, + "nr": nr_unit, + "nt": nt_unit, + "nz": nz_unit, + "spacing": 1, + "T0": 300, + "tube_k": "rigid", + "eta": 0.1, + "tube_mult": 58.5, + "ass_tube_per_panel": 2, +} + + +def calc_fd(Re, tube_eta, tube_Dh): + """ + Calculate Darcy friction factor from: + Zigrange and Sylvester 1985, A review of explicit friction factor equations. + J of Energy Resources Technology) + Equation 13 + + Args: + Re (double): Reynolds number + tube_eta (double): tube roughness + tube_Dh (double): tube hydraulic diameter + Returns: + fd (double): darcy friciton factor + """ + if Re < 4000: + # laminar flow + fd = 64 / Re + else: + fd = ( + 1 + / ( + -2 + * np.log10( + (tube_eta / tube_Dh) / 3.7 + - 5.02 / Re * np.log10((tube_eta / tube_Dh) / 3.7 + 13 / Re) + ) + ) + ** 2 + ) + return fd + + +class InterfaceUnitTests(unittest.TestCase): + """ + Run unit tests on some funcitons of interface.py + """ + + def test_file_not_found(self): + flux_filename = "no_file_existing.csv" + r = 10 + c = 17 + data_shape = [r, c] + z_offset_in = 0.0 + z_data, theta_data, flux_data = interface.read_month_day_hour_flux_file( + flux_filename, data_shape, z_offset_in + ) + print(theta_data) + self.assertTrue(np.array_equal(z_data, np.zeros([r, 1]))) + self.assertTrue(np.array_equal(theta_data, np.zeros([1, c]))) + self.assertTrue(np.array_equal(flux_data, np.zeros([r, c]))) + + def test_fake_file(self): + # Create a mock CSV file content to test the parsing + with tempfile.NamedTemporaryFile(delete=True) as f: + f.write(header) + f.write(fake_data) + f.seek(0) + z_data, theta_data, flux_data = interface.read_month_day_hour_flux_file( + f.name, [4, 8], 0 + ) + # spot check data + print(z_data) + self.assertEqual(z_data[2], 2000) + self.assertEqual(theta_data[-1], -180) + self.assertEqual(flux_data[2, 5], 1) + + def test_apply_tube_flux_bcs(self): + # test applying BCs via interpolator + times = np.array([0, 1]) + tube = interface.create_tube(tube_dict_unit, times, []) + tube_absorbance = 1.0 + tube_rec_theta = 12 + z_pts = np.linspace(10, 20, nz_unit) + theta_pts = np.linspace(15, -15, nt_unit) + data = np.array([[(2 * ix + 4 * iy) for ix in z_pts] for iy in theta_pts]) + interpFun = RegularGridInterpolator( + (z_pts, theta_pts), + data, + method="linear", + bounds_error=False, + fill_value=0.0, + ) + tube_rec_theta = 4 + interface.apply_tube_flux_bcs( + tube, tube_absorbance, tube_rec_theta, [interpFun] + ) + # spot check tube BCs + self.assertAlmostEqual(tube.outer_bc.data[0, 0, 2], interpFun([0.0, 3.0])) + self.assertEqual(tube.outer_bc.data[0, 3, 3], 0.0) + + def test_calc_fluid_velocity(self): + mfr = 10 + rho = 5 + Dh = 3 + vel = mfr / rho * (4.0 / (np.pi * Dh**2)) + self.assertAlmostEqual(interface.calc_fluid_velocity(mfr, rho, Dh), vel) + + def test_calc_reynolds_number(self): + mfr = 28 + rho = 5 + mu = 5 + Dh = 3 + vel = mfr / rho * (4.0 / (np.pi * Dh**2)) + Re = rho * vel * Dh / mu + self.assertAlmostEqual(interface.calc_reynolds_number(rho, vel, mu, Dh), Re) + + def test_calc_fd(self): + Re = [2000, 10000] + eta = 0.2 + Dh = 0.012 + for iRe in Re: + fd = interface.calc_fd(iRe, eta, Dh) + fd_test = calc_fd(iRe, eta, Dh) + self.assertAlmostEqual(fd, fd_test, places=12) + + def test_calc_friction_p_loss(self): + mfr = 0.001 + rho = 50 + mu = 0.004 + delta_l = 1 + p_loss = interface.calc_friction_p_loss(rho, mu, mfr, tube_dict_unit, delta_l) + test_val = 26.09267992590423 + self.assertAlmostEqual(p_loss, test_val) + + def test_calc_manifold_bend_loss(self): + mfr = 0.001 + rho = 1550 + mu = 0.004 + manifold_tube_dict = { + "od": 21.3, + "t": 1.25, + "bend_radius": 18, + "bends_per_panel": 4, + "eta": 0.01, + } + num_bends = 4 + + p_loss = interface.calc_manifold_bend_loss( + mfr, manifold_tube_dict, rho, mu, num_bends + ) + test_val = 0.12351409567895541 + self.assertAlmostEqual(p_loss, test_val) + + def test_create_receiver(self): + num_days = 1 + times = np.array([0, 1, 2]) + period = 2 + panel_k = "disconnect" + num_panels = 12 + results = [] + rec = interface.create_receiver( + tube_dict_unit, num_days, times, period, panel_k, num_panels, results + ) + self.assertEqual(num_panels, rec.npanels) + for panel in rec.panels.values(): + self.assertEqual(tube_dict_unit["ass_tube_per_panel"], panel.ntubes) + + def test_set_tube_pressure_bcs(self): + times = np.array([0, 1]) + cyclic_times = np.array([0, 1, 2, 3]) + tube_pressure = 12 + pressure = np.ones_like(cyclic_times) + tube = interface.create_tube(tube_dict_unit, times, []) + tubes_dict = {"0": tube} + interface.set_tube_pressure_bcs( + tubes_dict, tube_pressure, pressure, cyclic_times + ) + self.assertTrue(np.array_equal(tube_pressure * pressure, tube.pressure_bc.data)) + + def test_set_and_downsample_tube_temp_bcs(self): + times = np.array([0, 1, 2]) + cyclic_times = np.array([0, 1, 2]) + num_cycles = 1 + tube = interface.create_tube(tube_dict_unit, times, []) + tubes_dict = {"0": tube} + data = 245 * np.ones([tube.ntime, tube.nr, tube.nt, tube.nz]) + tube.add_results("temperature", data) + # Test 1: inlet T + interface.set_and_downsample_tube_temp_bcs( + tubes_dict, + True, + 245, + cyclic_times, + num_cycles, + "3d", + "", + ) + self.assertEqual(245, tube.T0) + analysis_type = "2d" + locs = ["max_T", "max_avg_T"] + for i_loc, loc in enumerate(locs): + tube = interface.create_tube(tube_dict_unit, times, []) + tubes_dict = {"0": tube} + data = (i_loc + 245) * np.ones([tube.ntime, tube.nr, tube.nt, tube.nz]) + tube.add_results("temperature", data) + interface.set_and_downsample_tube_temp_bcs( + tubes_dict, + False, + 245, + cyclic_times, + num_cycles, + analysis_type, + loc, + ) + self.assertTrue( + np.array_equal(tube.results["temperature"], data[:, :, :, 0]) + ) + + +class InterfaceRegressionTest(unittest.TestCase): + def test_simple_model(self): + multiprocess.set_start_method("spawn", force=True) + # Receiver filename to write to + rec_filename = "./test/Rec_Test" + # file extension used by srlife + hdf5_ext = ".hdf5" + + # GEOMETRY + # rec geom + num_panels = 2 # number of panels in receiver + rec_diam = 10000 # mm + rec_height = 12000 # mm + + # rec tube geom + tube_od = 21.3 # mm (Custom, Schedule 5) + tube_t = 1.24 # mm (Custom, Schedule 5) + tube_h = 12000 # mm (also rec height) + tube_spacing = 1 # mm + tube_eta = 0.01e-3 # tube roughness + tube_absorbance = 0.98 # absorbance of tube + z_offset = 0.25 # z offset for flux data + + # manifold tube geom + manifold_tube_od = 502.15 # mm + manifold_tube_t = 45.24 # mm + manifold_bend_radius = 2 * manifold_tube_od # used in p loss calc + manifold_tube_eta = tube_eta # roughness + manifold_bends_per_panel = 4 + # END GEOMETRY + + # ASSUMPTIONS + ass_tube_per_panel = 2 # number of tubes analyzed at each panel + outlet_p = 0.5 # MPa + + # tube assumption work + act_tube_per_panel = int( + np.pi * rec_diam / (num_panels * (tube_od + tube_spacing)) + ) + tube_multiplier = act_tube_per_panel / ass_tube_per_panel + print(f"Actual tubes per panel = {act_tube_per_panel}") + print(f"Tube Multiplier = {tube_multiplier}") + + # Structure assumptions + panel_k = "disconnect" # spring stiffness panel connect + tube_k = "rigid" # stiffness of tube connect + set_init_T_to_inlet_T = True + loc = "max_avg_T" # "max_T" + analysis_type = "2d" # "3d" + is_single_panel_analysis = True + single_panel_analysis_id = "0" + # END ASSUMPTIONS + + # BC/IC INFO + # flux data input and shape + flux_data_dir = "./test" + flux_data_shape = [25, 25] + tube_initial_temp = 300.0 # K + # END BC/IC INFO + + # ANALYSIS TIMES and CYCLE INFO + num_days = 1 + start_time = 6 + end_time = 7 + time_step = 1 # hour, can do nonInt + # stimes is index of hours to analyze, starting with first hour of flux BCs + stimes = np.arange(1, (end_time + 1) - start_time + time_step, time_step) + period = len(stimes) + # array of sample times of across number of days + times = np.tile(stimes, num_days) + # add time zero for initial conditions + times = np.append(np.array([0]), times) + # month and date from flux data + month = 6 + day = 20 + print("Times analyzed: \n", times) + num_cycles = 2 + # END ANALYSIS TIMES AND CYCLE INFO + + # TUBE DISC INFO + nz = int(interface.convert_mm_to_m(rec_height) * 0.5 + 1) + nt = 8 + nr = 2 + # END TUBE DISC INFO + + # SOLVERS # + # setup solver parameters + num_threads = 1 + # NOTE: using large tolerances to assure quick convergence + rtol = 1.0e-2 + atol = 1.0e-3 + params = interface.sample_parameters(num_threads, False, rtol, atol) + # thermal + thermal_solver = thermal.ThermohydraulicsThermalSolver(params["thermal"]) + # Structural solver + structural_solver = structural.PythonTubeSolver(params["structural"]) + # receiver system solver + system_solver = system.SpringSystemSolver(params["system"]) + # material damage model + damage_model = damage.TimeFractionInteractionDamage(params["damage"]) + # END SOLVERS # + + # MATERIALS + # fluid model + mat_fluid = library.load_thermal_fluid("32MgCl2-68KCl", "base") + # material models + st_mat_model = "elastic_creep" # "elastic_model" "base" + mat_thermal, mat_deformation, mat_damage = library.load_material( + "740H", "base", st_mat_model, "base" + ) + # END MATERIALS + + # structure output + save_struct_files_to_vtu = False + st_filename = "structure_" + st_mat_model + tube_filename = "tube_" + st_mat_model + + # FLOW SYSTEM SETUP + T_out_target = 720 # target outlet temperature, unit: C + # tolerance for mass flow optimization + # NOTE: making very large for test for single iter convergence + pct_err_outlet_temp = 100.25 + panel_flow_path = [["1"], ["0"]] + mass_flow_per_path = np.array([650, 650]) # kg/s + T_in_per_path = np.array([500, 500]) # Celcius + use_cycle_reset_heuristic = False + save_heat_to_vtu = False + # END FLOW SYSTEM SETUP + + # START ANALYSIS WORK + # rec tube data structure + tube_dict = { + "od": tube_od, + "t": tube_t, + "h": tube_h, + "nr": nr, + "nt": nt, + "nz": nz, + "spacing": tube_spacing, + "T0": tube_initial_temp, + "tube_k": tube_k, + "eta": tube_eta, + "tube_mult": tube_multiplier, + "ass_tube_per_panel": ass_tube_per_panel, + } + # manifold tube data structure + manifold_tube_dict = { + "od": manifold_tube_od, + "t": manifold_tube_t, + "bend_radius": manifold_bend_radius, + "bends_per_panel": manifold_bends_per_panel, + "eta": manifold_tube_eta, + } + + # STEP 1: create a receiver + rec = interface.create_receiver( + tube_dict, num_days, times, period, panel_k, num_panels, results=[] + ) + + # STEP 2: read and assign flux BCs to receiver tubes + # read flux files, make interpolator functions by hour and save + # NOTE: not doing multiple days right now + # NOTE: can we assume zData and thetaData are constant in time??? + flux_interpolators_by_hour = interface.get_flux_interpolators_from_data_files( + times, + start_time, + end_time, + month, + day, + flux_data_dir, + flux_data_shape, + z_offset, + ) + # Take flux interpolators and write flux BCs for each tube + # flux is sampled at tube centerline and applied smoothly + # across sunSideSurface via cosine + interface.calc_and_write_tube_flux_bcs( + rec, + ass_tube_per_panel, + num_panels, + flux_interpolators_by_hour, + tube_absorbance, + ) + + rec.save(rec_filename + hdf5_ext) + + # STEP 3: add flow paths to receiver + # add flow paths to receiver + interface.set_rec_flow_paths( + rec, panel_flow_path, mass_flow_per_path, T_in_per_path + ) + rec.save(rec_filename + hdf5_ext) + + # make solver + solver = managers.SolutionManager( + rec, + thermal_solver, + mat_thermal, + mat_fluid, + structural_solver, + mat_deformation, + mat_damage, + system_solver, + damage_model, + pset=params, + ) + + # STEP 4: optimize mass flow rate to match target outlet temp + # now optimize mass flow rate for each path + # this takes a long time, can skip if it has already been done + run_mass_flow_opt = True + # Heuristics + if use_cycle_reset_heuristic: + solver.add_heuristic(managers.CycleResetHeuristic()) + if run_mass_flow_opt: + interface.optimize_mass_flow_rate_per_path( + rec, + rec_filename, + T_out_target, + pct_err_outlet_temp, + solver, + save_heat_to_vtu, + ) + # Saving optimized results to base filename + rec.save(rec_filename + hdf5_ext) + else: + rec_filename = "./Rec_Test_ht_iter_2" + rec = receiver.Receiver.load(rec_filename + hdf5_ext) + + mass_flow = np.array([path["mass_flow"] for path in rec.flowpaths.values()]) + print(f"Mass flow: {mass_flow/3600}") + + # STEP 5: calc pressure loss from flow + # calculate pressure loss from + # pipe friction + # head loss + # manifold loss + flow_path_p_loss = interface.calc_p_loss_from_flows_temps( + rec, tube_dict, manifold_tube_dict, mat_fluid, outlet_p + ) + inlet_p = outlet_p + interface.convert_Pa_to_MPa(flow_path_p_loss) + + # STEP 6: solve structure and life receiver + struct_output_dict = { + "save_to_vtu": save_struct_files_to_vtu, + "st_filename": st_filename, + "tube_filename": tube_filename, + } + interface.run_struct_analysis( + rec_filename, + set_init_T_to_inlet_T, + loc, + analysis_type, + is_single_panel_analysis, + single_panel_analysis_id, + inlet_p, + outlet_p, + num_cycles, + solver, + struct_output_dict, + ) + # END ANALYSIS WORK + # do some checks on hdf5 files + file_gold = h5py.File("./test/Rec_Gold.hdf5", "r+") + file_test = h5py.File(rec_filename + hdf5_ext, "r+") + # flux bcs for a panel/tube + data_name = "/panels/1/tubes/0/outer_bc/data" + print("Flux BC") + self.assertTrue( + np.allclose( + file_test[data_name], file_gold[data_name], rtol=rtol, atol=atol + ) + ) + # fluid temps for two tubes + print("Fluid Temp") + data_name = "/panels/1/tubes/1/axial_results/fluid_temperature" + self.assertTrue( + np.allclose( + file_test[data_name], file_gold[data_name], rtol=rtol, atol=atol + ) + ) + # temps for a tube + print("Temp") + data_name = "/panels/1/tubes/0/results/temperature" + self.assertTrue( + np.allclose( + file_test[data_name], file_gold[data_name], rtol=rtol, atol=atol + ) + ) + # NOTE: in interface.run_struct_analysis we make this a + # single panel model, and don;t save it in this run + # so structure results are not available + # all else matches and run_struct_analysis finished, + # so call it cleared From 54373780a33995a6662980313fac5fd4c2725747 Mon Sep 17 00:00:00 2001 From: Bryce Mazurowski Date: Thu, 10 Jul 2025 11:32:35 -0500 Subject: [PATCH 5/7] fix: small fixes and updates to interface.py --- srlife/interface.py | 128 ++++++++++++++++++++++++++--------------- test/Rec_Gold.hdf5 | Bin 115736 -> 115736 bytes test/test_interface.py | 24 +++++--- 3 files changed, 98 insertions(+), 54 deletions(-) diff --git a/srlife/interface.py b/srlife/interface.py index 0a02b4c..3ac0a95 100644 --- a/srlife/interface.py +++ b/srlife/interface.py @@ -178,7 +178,7 @@ def apply_tube_flux_bcs( list of interpolating functions for heat flux on tubes """ num_steps = len(flux_interpolators_by_hour) - tube_theta = np.linspace(0, 2 * np.pi, tube.nt + 1)[:-1] + tube_theta = np.linspace(0, 360, tube.nt + 1)[:-1] tube_z = np.linspace(0, tube.h, tube.nz) # add extra time step for time zero tube_flux = np.zeros([num_steps + 1, len(tube_theta), len(tube_z)]) @@ -194,12 +194,12 @@ def apply_tube_flux_bcs( # apply across crown so only sunSide gets flux for i_theta, theta in enumerate(tube_theta): # for each thetaPt around tube - if theta > 270 or theta < 90: - tube_flux[i_hour, i_theta, i_z] = rec_flux[0] * np.cos( - theta * np.pi / 180 + if 0 < theta < 180: + tube_flux[i_hour + 1, i_theta, i_z] = rec_flux[0] * np.cos( + (theta - 90) * np.pi / 180 ) else: - tube_flux[i_hour, i_theta, i_z] = 0.0 + tube_flux[i_hour + 1, i_theta, i_z] = 0.0 # make bc object for tube with heat flux data # NOTE: this is done here and not in cerate_tube because if we initialize # bc data to np.zeros, flux will always be zero unless we reload from a file @@ -326,12 +326,10 @@ def set_rec_flow_paths(rec, panel_flow_path, mass_flow_per_path, T_in_per_path): for tube in panel.tubes.values(): # for each tube in panel times = tube.times - data_shape = tube.outer_bc.data[:, 0, 0].shape # initialize mass flow for path mass_flow = mass_flow_per_path[i_path] # convert mass flow from kg/s to kg/hr mass_flow *= 3600 - mass_flow *= np.ones(data_shape) # set inlet temp T_in_flow_path = T_in_per_path[i_path] * np.ones_like(times) rec.add_flowpath(flow_path, times, mass_flow, T_in_flow_path) @@ -416,6 +414,9 @@ def optimize_mass_flow_rate_per_path( check_outlet_temps_and_step_mass_flow( rec, path, path_key, T_out_target, pct_err_outlet_temp, breaker ) + for panel_key, panel in rec.panels.items(): + for tube_key, tube in panel.tubes.items(): + tube.write_vtk(f"ht-tube-{panel_key}-{tube_key}") if all(breaker): print(f"Converged!!! solved in {i_opt} iterations") @@ -674,7 +675,35 @@ def calc_p_loss_from_flows_temps(rec, tube_dict, manifold_dict, fluid, outlet_p) return flow_path_p_loss -def set_tube_pressure_bcs(tubes_dict, tube_pressure, pressure, cyclic_times): +def update_tube_pressure_bcs(rec, inlet_p_per_path, outlet_p): + """ + Set tube pressure BCs based on pressure loss calc after thm solve + + Args: + rec (Receiver): receiver object to set the pressure for + inlet_p_per_path (np.array): array of inlet pressures for each flowpath + outlet_p (double): outlet pressure, same for all paths + """ + # set tube pressures and temperatures + for path_key, path in rec.flowpaths.items(): + # for each flowpath in the model + tube_pressures = np.linspace( + inlet_p_per_path[int(path_key)], outlet_p, len(path["panels"]) + 1 + )[:-1] + for i_panel, panel_key in enumerate(path["panels"]): + # for each panel in flowpath + for tube in rec.panels[panel_key].tubes.values(): + # for each tube in panel + # get tube times + times = tube.times + pressure = np.ones_like(times) + pressure[0] = 0.0 + tube_pressure = tube_pressures[i_panel] * pressure + tube_pressure_bc = receiver.PressureBC(times, tube_pressure) + tube.set_pressure_bc(tube_pressure_bc) + + +def cycle_tube_pressure_bcs(tubes_dict, num_cycles, cyclic_times): """ Set pressure bcs on internal face of tube from flow @@ -682,12 +711,25 @@ def set_tube_pressure_bcs(tubes_dict, tube_pressure, pressure, cyclic_times): tubes_dict (dict): rec.panels[panel].tubes dictionary tubes from a panel tube_pressure (list[double]): pressure values for tubes in a given panel with height pressure (list[double]): properly sized array of pressure data to define bc - cyclic_times (list[int]): list of times across all cycles of analysis + times (list[int]): list of analysis times """ for tube in tubes_dict.values(): + # NOTE: as soon as I cycle results + # ghost temp, fluid temp, fluid velocity results are no longer good + # so lets delete themo + try: + tube.quadrature_results.pop("ghost_temperature") + tube.axial_results.pop("fluid_temperature") + tube.axial_results.pop("fluid_velocity") + except KeyError: + # We do not need to do anything here, results are already gone + pass # for each tube in panel - tube_pressure *= pressure + press_bc = tube.pressure_bc + pressure = press_bc.data + tube_pressure = np.tile(pressure[1:], num_cycles) + tube_pressure = np.append(0, tube_pressure) tube_pressure_bc = receiver.PressureBC(cyclic_times, tube_pressure) tube.set_pressure_bc(tube_pressure_bc) @@ -722,13 +764,13 @@ def set_and_downsample_tube_temp_bcs( T = tube.results["temperature"] if set_init_T_to_inlet_T: T[0] = inlet_T - T[-1] = inlet_T tube.T0 = inlet_T - _, r2, r3, r4 = np.shape(T) - T = np.reshape( - np.append(np.tile(T[1:, ...].flatten(), num_cycles), T[0, ...].flatten()), - [-1, r2, r3, r4], - ) + _, _, _, r4 = np.shape(T) + # this will tile stack all temp results num_cycles times + # It does not account for cycle heuristic + T_0 = np.array([T[0]]) + T = np.tile(T[1:], (num_cycles, 1, 1, 1)) + T = np.append(T_0, T, axis=0) tube.results["temperature"] = T tube.set_times(cyclic_times) T_3d = tube.results["temperature"] @@ -777,6 +819,7 @@ def process_single_panel_analysis(rec, struct_output_dict, single_panel_analysis # update model to use just one panel single_panel_model = rec.panels[single_panel_analysis_id] rec.panels.clear() + rec.flowpaths.clear() rec.add_panel(single_panel_model) @@ -803,8 +846,6 @@ def run_struct_analysis( analysis_type, is_single_panel_analysis, single_panel_analysis_id, - inlet_p, - outlet_p, num_cycles, solver, struct_output_dict, @@ -825,15 +866,12 @@ def run_struct_analysis( full receiver single_panel_analysis_id (string): name of single panel to analyze if is_single_panel_analysis == True - inlet_p (list[double]): list of inlet pressures for each flowpath. This - list will be updated based on results from this function (Pa) - outlet_p (double): pressure at flowpath outlets (Pa) num_cycles (double): number of times to repeat load cycles for the analysis solver (managers.SolutionManager): system solver for the receiver struct_output_dict (dict): dictionary with info about output for st files save_to_vtu (bool): save files to vtu output - st_fname (string): hdf5 filename for structural output - tube_fname (string): filename for vtk output of tube analysis + st_filename (string): hdf5 filename for structural output + tube_filename (string): filename for vtk output of tube analysis ToDos: rename loc to better name, make it and analysis_type bools? @@ -843,32 +881,29 @@ def run_struct_analysis( rec_struct = receiver.Receiver.load(rec_filename + ".hdf5") rec_struct.days *= num_cycles inlet_T = rec_struct.flowpaths["0"]["inlet_temp"][0] - # times we actually analyze receiver for thm - analysis_times = rec_struct.panels["0"].tubes["0"].times[1:] - pressure = np.ones(len(analysis_times)) - # create cyclic time for life - cyclic_times = np.tile(analysis_times, num_cycles) - pressure = np.tile(pressure, num_cycles) - for i_cyc in range(num_cycles): - cyclic_times[ - i_cyc * len(analysis_times) : (i_cyc + 1) * len(analysis_times) - ] += (i_cyc * rec_struct.period) - # initial steps for pressure and time - cyclic_times = np.append([0], cyclic_times) - pressure = np.append(0, pressure) # set tube pressures and temperatures - for path_key, path in rec_struct.flowpaths.items(): + for path in rec_struct.flowpaths.values(): # for each flowpath in the model - tube_pressures = np.linspace( - inlet_p[int(path_key)], outlet_p, len(path["panels"]) + 1 - )[:-1] - for i_panel, panel_key in enumerate(path["panels"]): + for panel_key in path["panels"]: # for each panel in flowpath - set_tube_pressure_bcs( + # times we actually analyze receiver for thm + first_tube = next(iter(rec_struct.panels[panel_key].tubes.values())) + analysis_times = first_tube.times[1:] + # create cyclic time for life + cyclic_times = np.tile(analysis_times, num_cycles) + for i_cyc in range(num_cycles): + cyclic_times[ + i_cyc * len(analysis_times) : (i_cyc + 1) * len(analysis_times) + ] += (i_cyc * rec_struct.period) + # initial steps for pressure and time + cyclic_times = np.append([0], cyclic_times) + pressure = np.ones_like(cyclic_times) + pressure[0] = 0.0 + + cycle_tube_pressure_bcs( rec_struct.panels[panel_key].tubes, - tube_pressures[i_panel], - pressure, + num_cycles, cyclic_times, ) set_and_downsample_tube_temp_bcs( @@ -890,8 +925,11 @@ def run_struct_analysis( solver.solve_structural() if struct_output_dict["save_to_vtu"]: save_structural_results( - rec_struct, struct_output_dict["st_fname"], struct_output_dict["tube_fname"] + rec_struct, + struct_output_dict["st_filename"], + struct_output_dict["tube_filename"], ) + rec_struct.save(rec_filename + "_struct.hdf5") def create_receiver(tube_dict, num_days, times, period, panel_k, num_panels, results): diff --git a/test/Rec_Gold.hdf5 b/test/Rec_Gold.hdf5 index 1fb9fc4f1af9d0ef7b9ef4e7c04e2200a4311692..23ef1c70bdf2643049f612c95e86350fb7cb4051 100644 GIT binary patch literal 115736 zcmeHQc|c9sAAd6|^SO#}d9)#DB#;8R89HS)eRn+uqtvQrHo+HiPmY_!(R(yj{H92m2}%0|yWF z_jd79(U0ms^iM0S_VW{wS_wgJ$K%p?E4$6F#!1s(1ikaEY zK;3SD%Sg3|n%befka`5%@95~R+-52$%q^GBrR2sdk4U+=fs|Pr{rT$)z z{+GVPafy1G6=w1h7;m|a%f!m8e%4uo-P8U5@-4uROY|pa@E^1Q8!ZY&Be;QnNw9*J zrG2;V?er9d-mO;S*E5PpHwvx~ZeVx7^%WVg?qEIM-?bFuv#;LKruz?STfhKDyRZiW zg<^oGuiH>BFLyUTrYg>32B5i*r-vt<_7q@N+9T|*Njp};bOLjtLSgIf;@4sD@BnOI z3F>$796ZEksQbWyo^GD*UVb)1y+%-;o+I4hdvqtBei8T{`xois^@q|g3Up8|)C=vv z^DVA{uMyn;P`)LgK2N`-2KuG6q>uA5_A8FR(m)6A}FkQh|24?E_kT%_`!1H#E}~ z1;hgUJ|ONps9=8>&EPlO)m_@N#(f7{{6c1l`+)EK4jdWjo5Qs&w@>YTgE7QbpSRb?#vRoAY0x>+9+#%L!(5ne;!T+$bDM9hYO;FhXP>$Gb>eGj0e-&0 z_3J$N50b-uP2Bgy@lLRUmL;2S^xo}fdxPI1UEYtiBwt{v>(S_^WHanHeC`-uiR*^nIhmF z3PmCt7HA$tVtXcTYDZSuf%K_0*asXJ@x$5md4Ecs@m`g2zYn%|rukYTn-zx9X`Naw$5)ZDpSAd>!*&&{#cwyfgV zrBQ`4eXX@EOBQzN_^4uEYi&zOKIxCQGrOz%Ur;SjqWYC*p-ZEcj%_NuR$h!!MwB$m zuo|#RS$S@mEl%B5DbtGGNg2^HO3B;Df2&jf?O$(H_N?jC_Hg(H<(NU1rR{=Zm89Ws zhL=ahDg)*VZztB+kgc6uemsBOHM78skl87d6dRO*Ew=XUGJn0&)%eg3mp|7l?Yow8 zXy1OlQq?lStK;I>Z2h3Wq~1~2;_=G8)fQ~s}T2 z!^dyU)=nobrTXIab3NYmho*VF^4V)UD#l$?Vu^OwRG+_(`dFkDD(7J! zk_JmpKAnG6By}B^ojJB6h}16S6IW<{FsXUKdQP(uA;c`Obikk=r<41ZF580kc_@oT zkak15sg&QR9vsm9gi86<(U_iFI(sS$`&K%8d8?;#?8&0zPn92}oaw!C^E8)1N=@%~ z{BuWc`wj7yOjCS5LA`Zl!LR`CH#ma2_}ps;&#h<$KldT@Cjl;g!+qSC7BoW(X{qQ3R_FI&qd5D2 zAejPRP1C%{zqjQ!F9x(@rbNFh@P;j52H@vKGtgfJvjdvLLp9Xb&xxNMGv61v>ex(U zsmXu&eGxT*L?5CnsReX)(D{Kh5I&(FY-;S}^&Q=3J27@ttA{@e}48HwK@<6&@cJMdX5UJfU@1@rpU!vDPmhEYrCXj8) zUeTML2&4r`3cp9+AI7{oa^@FPk!%<_`nl9@;hwkirFIQte~Cz~D3UkRtlFHYDUzu_ zJ`|qT5lN{x^{S?}6v=m&s{N{R5Q*NdS_U4p-Jihk3h+w;esRFB4)7CzUp(;J0s5T- z{nA0dG|=x7=rk9m80Kc2S zFCF+@0sZQMetw{zKj?Q0^qUI$MSy;f+r`|J>|%d$!l{?E|LXaeY&=5y@BESriC!v! zbT4E7+G{k`r{gx8O|<_ijWz5c$@wb%yX#WBTi0ieUPt?{+lzM1mec<0Gy8#Z8tuP~ zU=!=9w4XX|O>^HOkbM#R72lko<(i6bcF=w28sD}WA+@V(czM$w z6+}|ZVtvX6Gm*5fcEIu!?JxaDhP#H*{@77wZJEi#MG~bytCqB1J)%c~D77nhxX+8t z?jreN$}^)BPm#D>e>wR}1*n0DcvKpDpk! z3jDTl{Q81^l|Vmh&~GW|R~Yo$&goaL@hhby=gj>(52tmee%o^4;_Ys9+`St1zxsg! zi7q{3na^yxe}DVT6#H0#Ol{zjF4?7Oi^r{`cJ1eVeQ{-bflN5w*x(`^7sY;9JjOOa zAUk^8*|3t%XXOJL1h zu1Z+HK3Ap+r0Xo_Xs<|tc)7d&sk%rc!CysN?b=G~6)w8z&0&#T3RqUr<+(`O6mxGawR-^ke1YFjz%Lp26*|KB zH3xnVfL}1^w*&O63i@>b{f2>l6*>LBI$~ih$!QeV|H;ZO0{Oi}#uF7CKlKh4yOd1( zulSzXHBILTWX6Iyu0=NpB>dQ=2a;Wm&brcCYB$rNyu)4EfAig)*&Px3MY2V*fHpo|3YGZ-i-pf^;7&S9(i0bNoF&5z8%OUJe z6Q{DOIB@WA&jC5Sm$UmL+y@VJ^Yj~;`JQf`wRyIy?DHLKp~uoO-sNxb`;(d1a_2Ke z9VSc^CsP1-{b?`xd@)DP{=771aKn5E|4qTZ)>9Nw?^BzA88Cm$a2UH`3mfgc-(&p@ z$ZJ6>`gsmf{~lq0i=T^h51K0oR>;fp6TWY&#n$m--KevJ&JNOm(D{MR4|IN@8wa|1 zKsOKQ=7G;|9+3JhkN(OfyLmRG`+q5Y0tXDH`}LyD>X79$-g5F~@z)JU(tY0R!=_gZ z5y+YA%Qn7VLic~CIV4N*A)fubF)b%ZbXGPB^qcv+K(0?~*~kB|KuXs4x~e!UkP4^A zM7iD;Nc9lCm;2rdWJB)|6cRm*)(eZD3;Zozf0-4xn%k3E%H2&exeO@O$kp%5I zeP(c3krZ*S*2cTONG81rc%2iSMIg^6kmnA_V+8UTf;@FWo;VHiR0aFo0Q(4FpA@i9 z0NCdO*yk4w_Ss-^%UROPypkO;x2Ew62LpS!(#!F6{m)!f-P9R5a^iH`?_v>F~*d2H}LLhfG-rrw> z#!F`CZ#$6_ofwd(7RYlKzamWwjs#t0pU*+Sha^tOh z+*NM;m5;;9jmPqFS-J69K29q)UdzXA<;HLMIIdj(!MHB(N80)iw#QO#d=lGZoEb}P z-aKN#oc{YYE?!S;Dw z%JpYn4)o`y4KDrBf#$XJYCGG|hUT?QKkDHfK=WFDtl;D1LGxPd%i7()-77Pnf~SN2 z?6|CJd2x_Psuv2Xo$M=;QI<(g&1f9|P9wu>+oy`e^~ZW~#w#-OD|kBS&mBOXz93H! z$fE>#`hz^^$2=Xm{>;;n>(4wLx&F-4k?YSqo&0v~ZNJV|il2Wy>tem;G~Vvh@@%oX zqXja#`^p9DX&z8uqp};vt)jrTISC#S%csd3XANfmkjLX}cvaTl* zcbhhC<~I?^nl+A2d%KF{X5{=kf75)SH`Xisn~Ngpro2!uCpvpU9v_frHOO-f z#N;B68`!5JXP=2+pK)Lxvs~=+lV2Zop29@4p1#dG3FOJn&V6dr-=B0x=vg_KB3Iz}v^P)uR@YeRlL-=U;^8JM2n`cJ4s)RqhU$?7o@i z8F>}4oKT+TJ2=~#&R@Qoz6Uuqc$5>J{u<;-19=3HM*w-=xaBMlFQ43Nx%o865A&k1 z)HXkg&y$jyFU99g$<3eQ^Qh$JQ=vSVSB0gv`Bi8i%(KE$+k7iN?@Dg|6`zMC*MBfC zi}xdK{Ri7)DL20e+hd*>mfGfv@p)r%^T+r+GP(I=d|sK{{4zezOm4mz+6VK_u+%pH zjL$=pn~%olrOEXl%v0n2NL&BG_IX{(^=DoV^yfYf86oWkh~!M+sgoD0M3S&+Yqaqs zy07N%b=rQ8NUSPVySy+lGar(tgZ?}&!80xX1kIzgFSWA%F_8!skLMJ@sV z=+5n2E4AN)^Gi9~{CQ;BpJ&E0H>LB?B)mBl7=}Oh3hkirJhFfJT(Ru)B+IN3?aqYphK6hGvv67xQdL8>S`h)#o)IA7YKgJ`H%&~9hngYP*Bj0OM7{qlIdpXMrv)_TKg0QdC^gIP@2MC_ z@%N;AxPB0qqqjmHH&)LjAZ)IH% zah}3bpPox8^?~;1utYUvUB%|b%x)7hKjwo^{CXHa`u#%T{N+r4d!qZxw7{qHtJKfJ zuX0@v$1~fZKV@1{_o*OQK+}3?I4-O9_`aG~zX!@~oZ@q9W~uFY#wKh*`>f70&ctnw zpODr5{1uK@XRe?4dIjxTd%u$VTw$TXzw*~}1wMC1nlVqnbN4f!7xmTqH$6ccP0v^O z{Hvxye`ry_L<b8A{0C?e--R~G!tJTV=cvC5 zclG{&dNA&tpyW2OIEV1#cPKtz(h*}O+&PIfgWL4QC-n~PJUDv7A zEoi$vT@zaEJT;y4tvu**7m7T{sM8>n9E#dD&&oWM)VHh{GmN%7 z{7dQ_sh#7gx9iWZm`Tzm^nYvrYbd!f$H}nrZCch}bL>>@Fw!F^XjWOhFcP?|Rz;tG z@Y@0WT!CLa@H-Fu{seyYfFJ5-1Nt=r{Z@j0D?q>Hpx+J9uX6nE^V0dMwo@nWEmcA! zqdkh(IlNyatB#nswGId(r=C^a|8t#SvVP6UDE%|RWPX<^qb0kzPl`JrwKFbv=jnp` z0!ehbuy)ikk$ije@@bP^LFAcnwSLMX!SwuDw~PK;gUL7hRt=KcnE=18z%K;&tp|Q* zfL|Q&TL}EZ4DRig_&HUs*(mvuNJ^G7O^Q4aM84B|rn=-GOb!j(dcIHf5OT^=wX0)H z2>B-ZZAXb;uUe)_Qopy0+q>6o`wSBNqw|F2#-YS+=D38F14GH8aur({Mun2wqpJDp z(f*kb8bqZ2n-@69Uh2P7-FocrJ0z4C9xwkgv2_?RP4KwYDJ+b{r8Qb$btH_u>~{CN zP-}09UrA2bc$>HiB^gK>Gub&@T1_qIGCU5?#|6MQ%*(dtcyAVt)OgErkLz5T8_1QB7C*Mygb!6d`c)b)9MFuAwA^%bez zGT_$&_+12kj=--x@CyfijeuVxqmHF5DZc~N`#&)37erhZgsJqJ2NP1^z>uRQLr7br z%h8G}A>>}=>k*-ArjyE77RPu}oR=9(92}d{b}o&g-8~o2Br(?GjP+ZDlC@q3!Ymhs zk_R5nyYHlhlB%y4->%R+j5NDc_=toLLE9)^-<-m}Y^+lLX4zXxr* zx;2axua(?=oqjm!U+Awh>pF##?*nJq{X*LX13zcr7YY2LIet}uUrpfmJ?NJJ`trv+Ig{T`>vp?2mDw)FgC=l_+hR4a&dUDWmIfVV-!e#pIxcC_Cf z_igFibrbc=1z)w3>{91m!eptP&yzVDrc|Tn0r>9!Y32)&B!_k%<#Q&8_)dO)?^bLu z34XEV^@xxV(z{q{1F4-G@GA%WE&x9p;MWcKr2)U+kl)i1-%9;gy+@b*55@$M$?Kd0 zVjP0W+EOlek{g7O&>3B)3a>)Qi$&gPiw{mGi>8^ck@$6u>J={a!_>1Ex|sJ3C2gO# zI@)VLohPH(SM{RvU%c_G(cOoJk;^MX=O@O6kuev#r%C7u03s))DLIUi>h3@iDXaxTRSS#Ja1Fw|4wu_2qx_k5{kaQ z5=_dC?c}{BK74S?7XC+qNb{ zo^1>!4+oT0dCv$T?P3PhklGo1%lIV$zumyE4DhoCe#L-aW#Fgjz1)!7H^}Y1Ja&=# z2~L;K%gj*V{`2qUW?%tZ*mpvM=)m9W@6dhs9z!p7xAhJp9i1Yh?o|&aqm~cbXml)?+`PXyc&%FqIih@WT#7I7 z>>C!ZaX_N;tGUUHM&!=?#BJ+WtbS>Q?+rGRA53pi_-QNsV4Z+L!PTB`@N;+3)e7 z@vEy|=W?Qh^2`Kzu7W)JAkRXO$Dj$5Cr*PrXrBjJ=kX~*yj=0=adHfta#dH zj}-seXSXV0Pa%q1=G>{9)pm+kC(b;$ZEz52r--;}VH!-P?YpwV`A9Ia9)Bt{Cp!8uw~bFfradY*J2i zmVrEnK%RLZPf?Jk1IRNSm5;N^jkoe~SGn<5J`RiW zXd92^X|}WV(ej6bNvXe*i;Y|uLY!)>aj^HAPR5@*6D`Hhc|9h!Id((JYw+m(pjAt+ z8RXKH@F5dyXOcFp=C*lVG?cXXd!V?wUnn_y_E~So^`XSld+rP=KY^#?VY;EOl;;rG zb#cQ14MWMeM5hGH*`ee_e3`_duR}?LpfWupx`vTI9@+#ypBqM^H{Z_4iOzVCXCTOP z2;^y@L7s0k$TJ=6^EKGV9_;fy*vA;`(+}+9rNKVF{nDpOb}p+NVcU)3JonhOwbQrM z-n(Kdj9f|WUHkF8_*zGUN%p$2YOS$<)*uDrs5HQtZ3^&f1X*QH#4=H)FM`}v!}*BSfYU$3Bj z```JM?NMgy{N-1<=J}mwQlsA$G=C}R`to^x_`76V_#HpmC$GpWkc?q zgUs_V5uN8$SuA5r&PX^ERc0#aG@nLL~?I)^R%AU zL8MOV^%KRO29bfqL$@Wa45q(FX|;3&jSnLqi?l-JJPbtAVCl)H^RJ4euH&*Z$9B+t zpHev8u1Nwq^O@Q-^+kysx#k zCE8V27EB9_E&XHq#eo-Abn%X+wRJxlUIAXl{rqZOrxsypFZYm9n+n1|x_hsY!3HX6py9^8RsK#5>BT zv2A`F-84uJQ+wcX9Lu4^q|U9*OD%@)?R?y#;o zz`EA&AM2W=gC_a&iXPgj3fK?J?C)vp$usZp6rABFEBHiWac60PX*0i3c9`;-~ZvfKeW%kg{Q!u zFDUnYx{@R)II&d7eV_jF@e~4JU1akj_OC)AWPZ#CpZN9A@bdeG!uc!q;`{3rv~RPo zhgqMqq`fEU4@>zt9kHKHZfK3usyo9aFPfAg}ksa~r3bI0{S- zeTBB?7h9ZUEkDoqA;arg?ayD~c#TcvzbV+)+WVE(^Re;kcL!cna~r&j`<@QWR6;jQIL_4jkx3F=~NBZg{ZLQA3hiS7U}YS ztR?vZQ(Y7OP}Djx{f$U^%#zm zqCLZL>RX%PIO+YuaGZJoj#H6t49Dpz;3|cD&2XH23>l8ooM?vQ)Tt64996j*Q_{Bo`B=z@62$V zy8p#+oP5eK9H;yC49BU?0fysb4meJWfL|kyUkl*p4g59$zd3;8bO3OUb^^}Pw}lvv zQypK1<20)o!*S{u!*HA&09R=m;3^#fT%~n@tJJ7G!*No*Oh8aFymdGaRSELm7_KaE|T(&e15qIr0RY zqeXyoRJ#hpaS8>TqhA46$p>(i_5!X_BH$_+s2Glu5^$9k1Fn)Q;5f|#94F&p496)D zaGcBl$0;0eoR$NQ(^TL$2>5LQe*1x62H-{^zewQM81(xC^cx8J`G9_NK|c%7F9h^+ z09++Az*R~CT%{3!tF#($l^XYBI8G)D8IDtZz;T)aI8OBd$LZ@ZhU1h3I8NsQ$7w&{ zI5`5pTEOou@EZpFHUmEk;8z0pHHGszYyjtI5#Su11RSRafOFJjF~e~xBQP8%d%#ti z0XR;*0avLF;3^#l9H%#c<8%#hm7W8R)9--etxX))j` z^#EL@>42-$5OAF00mtcQz;QYNI8I9e$Ele+!*SXKI8HBtUkdQM1^kr2Zx!(C4*VwN zg5#9L?HlCwULH_{?Ihs7!3R9`>Z*J}7U18@%|JKfuW*Nn-6%LXPEoM4_ut|;sW|(8 zK`!^Ns-kJXlHc#mbzc0A`M=q}Dy@!qezgEUFIs@0<}f>;+@A`^$;6WRqYlT3T3)wr zL)8AKvxCkKk|XH+K<5WKKhTW>-8`V12XynmhtC6&t@7ya8Q?hG1)QQn-VDdd)|ufr zU0K0!oYp;LI8KkO-;d*T0dSK%cQ72MYk-@S{F31~{btERn&;;n>(4wLx&F-4k?YSq zo&0t^18|Bu08Wt~;1sz5PLT)T6fFguqNf^gilzZ>l0V=k8Ky8Cryl_~=|{kEdIdO6 zi!|USc><1;2slo20LSSz;5fAd9H$q6<5X@-&NxnEK%QYB&m+KbasYX3Ie8XqkVgRf z)C2ptgM9{pecZr4(O{oP8tmf@xJi`(H>oe+Cfx(vq|tzz^fTZl?bLvq^bl~Iz5yJk zvLhLeQ!l`A@&O#D$AIH>L<5eK56DvuBm=}em zw)s(fo|N2tDL!vXZvGUXMrs5HQtZ3^&f1X*QH#4=H);;n>(4x$ zkG;vJo4x}OVgw{=Wc;W)|0co**bPhS(`4(G3ccJHs3sc*sflzlzS`kaOMO<2mu zak6X4ET6yqD!2V`p^I)kEEpprXDv0Yhi2S4-g&(rp4&M6cpNA1&-#5R>4*2v>+j}& zuAg%5<$`lpd%u#4<0NSCFCkwU8}VxRbM`Zz7v{C zlb*WuCto?tCl&_RpIG4Q4P|@;^9xP56tOGV_r@?yYB~qVJepOoje&|lAN|#QJ-hD+ z_cOvWIUlN`|396y?>xQ;?3Ie<@;EV|lrwGLl`Ksjc6pRNk4?sV# zf@{aY?2a8s9|`k6)-zegK3TzdB~UL3+F?DjaQ2D7xOawn2RVBDct&ex3}Efw(Q^lS z9;`G5dd|JF))U}5>L+#|{oBxj~*Y1X8$o$y=4j3B>8eG2@qE0&#v2Xm7GoATepbBoCzdA%}l+4UyW_ zAGY**<4Xb=XJ59bahgDml3vl9o(SYO>!k2|1|o6$`^cGJOhxi^^yueOyG486&X?MS z#r_hJT2UnX1*AhfEFx>9A)b2R&I|}@qfZtZ&*AVz60zboD zjNkd6+#WTM?Abl(&E5V2X*FK z6|S|uj3no}rgq(Cz7$CH!|ol17ZHg|#L@D-N{M7qA(xjP)kU(aTZp7?xanxA zU2{ch4XNE^i>o#j+KS{_CATKAKZs=N=I9g`JCO`r{P?7at4Myoe*JB!u_Cz`WgR57 zds~z7ivfP-z|R%EB(K+V#9KYxFwW zf2J?nHCs;mum79}%4xL!hD|rIo=*F3bFCpI*UOc_l1k17{oF-aipmUkJsW1&ENJ8VsJ-y)Fm5&IS2oS^-- zrsA6&G|#I3x2;Ad%;e!B zxl?FXEh#@MqDO)#wYzYr&x_6OBB?RunNf1PT01%iHtpr0Y=_kH76 zN=eQ^`*$8r>nsq*mJ1hecN0j*SHu2SKTsf5O3zs4Gg}~YZk?H8A1jbm^yO_6l z+)8TKVBXgkSGK2iJl@#gA{`geT^5hA4WN4TxU*p;ozJ@aHwbPWBM^Pni2G7I0{k3- zUpL^_6Zq8ve$l{hBk+sbdD>3mC-y&k_8c8|=R*7SxiVECX%WuRUXcR1HX#0a;cWu( zA;zy?ou=b_@2$Tje$SE{1xoEYPPtcQNgI)TU9m>!n)V_YI(KTt=T0ISb+}N*H9wIA zEtukHwooJu>n<>t`r&-Wk#MOU8P{*jb$Wlf{v(^Ji$qf9t7xlTTj~5+xag)ghegsR zU|B_%=OP(a)V;OT?ji8A27V^MuP^YM0{r}eUrXQ@3HnV3{jPw17eT+fpr1SFcLwyU zaKyq|lJn-4{!dnR5y{^Q@v&A<-Mup2u#1kmWSq^72*j*9}L~^2V^~6+>v8>gKYIua^j9vFMO2#fNzI z3&ymZAkit@IM8qA?*dsMwCv-5SRl0;cwJSTrSXe1W1?Jd3&dcC-phS&X?$S%l1Dkw z`PKPch(zan(+igx-<9Ora(hOGKpOwleO@O$kv!Ua`pn?6BDpcJS{v{BB1u$?e4P`W z*&vS|$g>3G`3B^f2=cs{&g6NjL7sA8pLSrM?qHwWV4pW&pB-SI1sd!VV{*$`(rZ>F zJ7R7vkhMP<*cIOVa`tBLl>2nQe)vnf15ZZ?p$2YOS$n$Y>#nfEVYd{^Koan@n=2`EjJ#`$ED@Q zr};Rw+;}w~x0V~f=Hu9MRDS{DE5;T=HpTDDd6aq^&fEzQf@-M`%{GoON|gZ}*cvaaRDK_a%-k~RVHoZlkdq)gC{rFNDk*wY6d%Y#~ z>$fK>BtP#WlGvS(t4etWJe{f&ANfmkE|jx5WnE7s&1~AVncqYt->-3W+S^qm10v_& z`J3ho^|fB%-&~}5D{U^+%Zbi@kf#I4;|cP3f;?M5o<|^0yasuOfqm|QeVTxM%7cBL zgMEgAeJ*RTkE>rFb)G^$v!1@qI??lyc6RPlo1UjM&T^#79BOaH*5A@^P<#J5^M_+= z*U@|(x3N#8JO$o9=B*yJknFR(_d5R~G~Z!-LbP)Sny(@bnC!lp<{4f5%5p+^n(wfr zqv`zRs|8ZZsllV1=%_%RZXnNAkS7S_DFX5&gFH<&$ivGg_gZc~4f4ahC@i(jkK*&B zpz&6#ru)A{)6qYl$&3K z?J-XbOKtPT_`ET>`D1(@ncRFbKCet}ei@%;;n>(4wLx&F-4`Pl31zt=Tall=Gb6e8hw zyiuGA0o?Vcy~Ohr@%O9oGv5EZG&2y-IPCrPKCKr2J-7ePuWVbfp2%N*m0O-);v6P5 z`rUuugE|R5r3LsrKm1*?GyIMpWy$OFqjhp-ANhIGU^>}0jlD^)v5pXuQw`t*7Rw6IDCV0%pl9sc0sX9((pIK%cEkI0rQ2o z6Kib9)=n-zc;((|3$|`@j8i5gn@&45CQkYE2g9U@s&UHRyNrb4{;Ic!vdQmzo9*}XP+t1ZuJf3Y9?GRA)r*;l z9?ErBsxMwY*W;aCtDL>Iqhj1OC6;IzT(uD&M7nQx(b}3X5&wfP)$QRSpfnPt<;T+f z2PxghFTeu1z4s@rWFCXhlV5yvXN&QiK9@83N80N z`pd^rFa+cs_Fjbjo8NI1OrmgQ%ar(^pK$(4{O$eq{@JD=e9rzpoOPcC@e{C=kE7tY zgjqg+`&F*%;Tqj~SP;gACLJ}chbr#e>Ac+#6r>YZ6bjTH`9#6F z^anREez~nbiY=@j#sTt`@pQuc>lV;k_LnA1e*C-SOD`;72Ixk4Ikx>>l zf0+L|Zq@m(mZtxrs`LJfZv#=_Kd7v?i)Hd$@_0T$z@$avfRf_y~m78CE{QFImhUd@F{AUa*k@ohG`&^!=!7ll0oP9cb zoc~Phk@Ncmzdq&mJ`uT#%^Nt63BGwuCD?sJvJXJL>31G|*YLjV8( literal 115736 zcmeHQc_0}vCg_6*PmBk@Q(*VYviwR-(T*GJesi;- zLwTSi=Es)(nwKe#K1tvmLffc3U^GbJ#dzH+k+hZgwu7Jf6MFEFX6p zFGui^`NQ$4iwpiuE62tw8P`%gO_sMAPl_kObGP$wnl*#RlLY~aVw7gTM+d^ytj38k zxEVHnOcm?G;d)E)j99)FJb7>gPcJ8XdpA2z+`jfAB1}!;)r&Bw78uW5;Aj`~y$JXZ zzr)uhZl{tcsxQr6Z}nZ5SshUMuX_!&Px$}W_W*fa;&`%V;|F^HT`fFbUv>wMOBxpt z(2g5pVkp9E@@2GIl!O>DZz8*$lY;iJq0g&g#~oabjCTR6@zWO{KEz}c;|ulz#(qIJ z1U%jhC(r4#+}!M@d!eS{eoTyQxrdX36X^DMtXshuL4W-Wk!8CRwol~ojO=W@hPinA zVtjkH{4ggM*CDg)?CqVVJK4E;4VmRO2OQIBjvf0w4kwbo8T&mxznC9y|6Ts(EDx?1 zZWs0gd0fdljw{9P|GUSPX3LZOGIj8`6p$bH%lN$F>ram5!Ta&`$?oQrlP5j;yU+*+7AfQG4jmVex+3sxIVqP+i1*C5@#YMpE zC-fOUZ}`4rQFr75S@wl^ef>aOFI4_Z|KDglC~PS70PY8pS$DXz)`I;Y!V|AG&P;t# zL0R^~t!rG$nUBgnUdR%&hP?0aRYG6jejtwB&gzQx%Ob;5TKaW4BX-+yIa{ES8GPO6 zEI|#oVB0Un(`7!vIbN9W;KA;fQ=-_um*#0p^0Ys*6OY*M7tVvfQ9Wo9&rVcd~( z>OQE&^$MywPooKwr{TyRzpgxuf18gDi~< zFWDGz2V|z;QKKJt=eoR^veSfaSpUDc%|y2dtX>qSAxF zs_IWQxDG6A|FY}r^aZuaW%t{T88cO1P-}9@WX2iZRDD4$tRMU%{fzeVcs?0cmIoA5 zDf4a~24_Q4sQJ#?a;E%biq_aJkv}(?@{I}nyrgRih09}mG`mT*F~Dh`x+P-?aFIF^ z14@IEDaGih5Fdb>UqSn~uLW22o3Z$Z&Z(4s+Yx^4*QQVrDNg${1CuGq$>!?A{gSEo z$_tw2E7d%nGS5`)V#ic!r^4j5A6BPO{T{T4J-8y7x_xy1vync@RNlNf>o0VviQD_A zWZW@@R7!DstJ%#~rcfhVI`Jx&CR1^Jd`p&lB~!Tp!PZYY*2F!Vl(6Nhd@5zz`|iL2 zfhm-F;JGGNi;}5cuUv;zx+hb{lGGd>g_^h(Epv?|WmBmePgiz}Seim@k}$aZd0sL# zd#hou7?)(~Mqiy4Nv&(*T5hmeY}+E058OTd-+>C~%H61m}>(y2*Dr2@SdWz;x-X+1Kt zF1E;|KA5gn>J4zl>hW?T0M4&xuU%3AXZ7&Im52T{ar<-TZ*SE+lRB;^Zo6tiI>qZ* zl(A!OI(4eQ>=3UG=~UUk6GyZb*2L*eySAo}cqUckpT8!5TskGb=-R4EpLEJ0)p6Ik zcInhH4gG0W^J?OPeJrnuF8CQ$MDw~RCHOi z=Cq;d)P_#u4^MMSr;_yYO+1>XQ!jLbWEEX%;_mkJl40 zor-bob&}sCow};}ct*~QnmBwt3GL720X#klI&R$IFAhmkcA!p3>7`@=+vUyu%a?(z zFq-qS5Li=@s#dRn)mC)qC#K>+?b>`|Tw7qVck1fL(f8$F%GFn%CjZqy#^C|{zQKZRF23(IW$#P|!yXc_p?|YD`5rcr@c-YWY4pYF z>KDo9w)*--Uqj?d9CwZ0v5m|DWWT7yw%?KM9k4qb9A|lYS&_YC_46W89K-6gb;U3F zUW6gY@}PwlctGd}VI05<;RnWn#X`p?eZYGxXa7wby)OjOeMfa3uicDKZ#uuy``~0g zoipXl60-z8E$VuH=v&EP`pAHHuU3HfWk%k949??|Y4ccsyWsQt@<{;qaj|Dm0l;ar z)=ZoOa6W@3j8>Eirb{ZDS3=xncOP^-lR=_u0M1vzc;`rfyWzDmEf3(fb*k#+1aN+x zPq&i;xHa+;3~tbH7Y6t6$ZLH^fLk^#@O3YM8xvqDp9FC65*?j~0o?69V*{UpdW!G= z$kg+=)f5KT*Y2pUCcrI9n=<*8IG@h{^I7CPfIGA3wU#`eEu*7 zcSv@{!;rs-`47yP*&gX?O0q#zdH zRzw2T!$MZkK;GWJ!xPxgfIsO1=+TG*OSb#ex zm1lGl;1o<$RL!>64)?4H`pFql5LuDKsP4;xpzGPv!>vPUEUT;||T5f1^*R>6Gvbb#x= zI>j;x;NIwj>SzMo(c2vuT=wvp46Z^nz9byr%metbR{(BjPu(yxfJ@u$^H(guE$#R6 zQWt=mdq{!71-`IkaMXRivjf0gQ!LrN1K@(@*rY21T*~AvAH4x?n#t0e6+?n)^AXYv z&Ue@-24`cl<9rW*ds(d7&lTV<+?g=xX$Xif^9~*L0IpBl&c*uxPWsJzU)a7o8~QW2 zRU0~wx&q>)`RoZje*?H1M|0fc0B)vvf%J>pe7Zbdp63W~FUlSV zuK#`Q_Lm0o>7tKnE5-m^Mx^a;n*mO;iTab)0GG7v3WMwMOp(FOe=w}f1mMJ|{$&XO zC+~jKsu#dj9ZW8s3vkkhEK{Bs1k?P}c?>SJ!#w6a_im5vB8S0w4E^vVA{W@vN1eAm z1vX-(zo|HQ4xhHaLS7kI&pDpbH^BbmJrJ(d(#5SIF0sd=c!={UZ?_lXvXcs}Ax`zc zxSkL{fhRF5Kx3aqHIa?E-O9rUnHN_hu@>wP=cP>m2T!g}A6@N(}D$h}jJ8 z;)*2#Gp0#DU~_4NYscbdS?8!j z+`Y-V?;);9f^vI^8(=>&2I8U@YVjcMHH*_2f{ypB=(Q5WwU-_D9OAT2XiGz!xcZva z5NEOCR0XW3&H!}0?z#xqA!CU&#A&QB%!jyrCNl3|1KZ>5!TAuExa~m!#EsQJI6YOi zeFtAuW`lB57TIM(+@Ac_MG)8K_g)SVcWYwvvk61+w#9bc# zAs6CCmAQ|DxRe1J`yeh(0pZfSA)I}${Cg0mee78(#0}94VQ?|6WEfm`-`u?rH?1SW z{oUtxQ8+GxKm0I_f0v4%-|POK>+`$gNOZ)k-;l-sf9&(S;}CpYp%ugf^*xt6jzDZm zf{+T!{xTjQpUc^4caS*sbrbDfyq#v$@?4Jg&#`lvHQmX}ulhM%`D=KMt(yHELWrPr zj4$~cRxP`_IeCmH9ow3>z;S_zer^(O(z)K3hNVd7G)Q@ObDn3jHAT z1H2%FaUhHXVH^mr17SZP><5JXz|YqYz;pJU{yM76zNi7-+a+lyK1v7*qA$u^9T1rl zL?^B9sMWg@pKh(_D?iwoPxqbSWU}ELpWYp{?^ytN-%lQX@P>s9Pq*#uC*%O0`^B4E z>;*iB`$-lL2RyHeCpDV_cm`?cW*-1N;VqBX%9GP%0mGvcJ#X+Kz+-!D#-l{QbHDG7 zoPmJHIFP4d3V8M%&7P73cr-$@YUMFd@nU$EeAaBY5%6R#tI}HucvRnZ?xzHJ?%(e? zRvYki{AeN;3U~~rW!B1LecYbm5i7HkaR)ph=M*m30-l5mY7bs_3?Oj&ilY1x_O@qL*HZv(Z?&Y5`qDb)hc=at$=5A?1EKq&+_T$_p58= zS@Ij=`E96bw?e?PdcR4J0>ERN_H>OG;CVP~)$5gjXW;&#-1~sXE3$STG7h+IbDfXF z$HVig*a|woO3t%#oo^-QUAfM`lJl@!=VQrvS+4W5uJh64yfoMOX>y*L>wGmiZ_Rc7nw-bx zI-gC>YjcesJkL$Wk)ZK|@!#9xc+RrVd!05AjK6$?=+CZT98WElo_mZ>%T3eTw-w}D z+}9o4D#5!te}b$Bj^{I#?yvFzPwRkEx~ChsPg!ohJ`V8gQC>Id4B#0RD5}1$3Isejj!7OI&+{{P>l*_e@8?P)v4E#t(YVzL zfTyqF8NMOl@vH1pZUA_29Fsgao_!_`z426!PrHuNyg3~3sLhHGJf6j;b7NPnFK-b{ zhjy3TTIyAuk3sU_cwVVEY4Hlcvt(W`ml78~U6%1qZwBD`6dM&04tQL8X5|F|9vsIc z502-Le;!d71$gcy{ppqlc&;gLH|PU+=1Obaw*owdJN8XB13Y9LlY9+jTV|RIlaFBA zyUT%ia36j1-lv=J_lj<%VkfTyn-{R-^b26;MBFrN1?ItT^^Glt`^oyhvO_%MrDx_q zo+lwgiXqRsE0&jG`Jp#XUx7SBj~+LJJW&H`=Ml3tB0) z6$vxgv$S+vst$n!?McAkSKh{wUPT`c4ov_&)v z@<=B?-vN2VPRi|uJby`NsX(4%-D~HuGeSJ==1*S@dFEf3lLC3dPMNKTJQ3=LVj<5g z#jp;LC#rMpJd=zOPf}YYhR6T99>cRDsWZdV({dHVqw~@M@hH~LL&gEuZLaw=_;@%k z3R^+*qsTlduK7}A-W1pTDKd|WYd#g3SH(5Iip;a(nr}tsU2)C7BJ;4g#t+WRBI8KV z_`!H=x#kyPJkArtR?vJgGH;A){ur4@#xDpOJZJT=UV$ zyfm)ygY(qLI1)5|F#da6uJKIPgKIpKJY3_Mo=G0A@l5jk z)9wE*pDD3H<}=locW4b4`2>Fw$O!*8QV-z#Di1qPu)>rlZ??D14D_>TE4+*GjE%38 zjm!V-n&s_f=V4_#9lOm5Xb}DXwmU$3f*)7oBL6gR3$K4lu76a$PL>$6zLl4sJM)9? z`qs5#;+_qDn6pi;M<)1sW+-zLu7?Ko<{Vkakn3JyKh(J%*}uH5Sj~PCW!`|8r$xgaiz<+K%^;(c>Mba?7!+f zV&tis`?=;9JF@Xdi+z0#@D=&>-hXSWzVX+cT?Y`i``@aEkhfuapw9efcKttGYn0z_ z&a(&ldvG5fKWJP)K#Rvwa=k&W<0zw}#|AOH;JA*XW|>`+A6l_fyz51?y76b0GLSEQuQERH1&c`m_aBe{2As$nh|9-nT1-`!8$o+Y{kmwFiD& zTowEjaaHelxCr?U$5XW@g`XOO2kIIRCH%i`J$_#OUcU#b?>fczt<_f0^^6tJfc8_Z zXI!n@oVoaG|C3wzdbK9b6IrgYZ3`Y(T<xi$0bzgz+JL8m4B;a-TtM3<^GVgE9 z+1Au`e}(UV85Zz6fx2`m3YMee{9T<2UJxpe;~8xU|2Gy7;PEG!ZLm9QZ z-9+Elr&FPD7nRd3KozJ&#uFD^$gz`r$ZG+Jv_0&B&Uv)emI*e_T%D zI{HGHydG=g=>USQpB6u_x_q&?XPaPrUhGq{1>+!@@Q7jidWbPuHu zHn-{d_DmQ(@tptrjsTaZJ!$13P|u*i>03krPA0;S!O3cV+78<{Az62xD!}>sPf>0I zaQ)w$9W{P-1f8JNEupl-I{JfOkYYEnXnIFzYX;}8?NntKHOpxs>iqW0^6-6pf}I=dQPae#ugMy$v33^{LiRHM)k<^NyVhy>dE? zUK?p9GG$f-turyq{kp%N2e`xK zLuEd?MbMI`B<{@xIQ{g#nd<=V#a~7YZhRk?aCkhsDI50;y&Xy`Uk^>rTo_L0Y*=}D zrg9`LqLQ{#>(V;+GbQrjR!A+RioWZ45ESS~zRwyk~eq{5C1>yAkuEyS`y&~zI8wakDySR=H7i;#| zGA)|!_*9m`HI*$K-Tgtko zM$_IYYhOaYJY94@6ykcfKJ=pEdMMqpWU^w_+;DoR!NKQEdqmQuPSdaMKDUlmn0+iI zKRKHAXsW{C4$jPHaG{R}4V(jT8XGhAdIQ{!cNZot?H);+(?d^c0Nnay^piUJq*(a~t4x>kU3-*Uv9rjZIxG3lO2F*ZuEs=v=OMxA`2I0ch5H4q= z%UX!j$#k6!ai!_o+#$~0b?F}vr`Y^$dx*=5K)Bb>P(9^u&+QCx)^-M;U_FQZULpo@ zhb>}qVLja|b@sw~s#zeMMgYQzP3!Im+qbzqb_B#3>CLr*?K`+qTpi+;DTKFyIQ~l+ zrhUe82-jqU!P=VuwQ6Q7U0TMGMvRgdC`#0lO6-Rb|J#8+Mka6S*m{f z2Ed6FDrHFn+|upN3M+`8)sC+eJ1264;Yz5NT~jvw*i z!Q-uVD!l@6B96OHL)^#_2zPAN8|HX7MklOS1v)$#&q}+_061~^V-+p{_q)AY!emfh zd7f6Jo)*xr+w6W1p1-D&RSeEm)!ji2;P^d7%OI|m>6P8C02jMxg5?x|^URAIJ`mvY zI?ZNqt7H+*?M5$4h$~ov5chti&O=zwR;rShVLe+Fpn7hcjBwAwY!xBS&H7g2ba=d#hW@a9BmTOo z4RM>D@AZT@Yj=bzy78K6U$bQrfAog!yHa_=7T}(ETJM?(aBW*Go-qmFcDkfZ9RP5J z<+GUf_4Y-$CGpMXr~q7x6N8&Voc>g!*UkW^bIBpq0^qt>oJslv;O>VZ+>{TX5*G@13KUBRCyiNIE#|{!6z(QMA0zgzU=*!q?^7#}X{;#2)2@lZ-A`lJ(rvF*I*S4x zAF21!ih!qi<$Mo6z|(!a)^RwhUE!X*4a^9Bf{4F_;%Qb%Rye=6>g2oTVW6O1Z665haGq!@x zH6|uz z=R{=8-gbaTJuqDV6{yd&PLlKV0grje=le#02gfnVgX1}6H?8gp+U2r(x}`a2SHRIo zjUj+%)zz@OL7-i;EmV&Lf_C9JCV6l?XTLQ2({FAVUHOTBMz?bWJ#%Bq9fLEGv{$R) z!~I>OXiwSm>h5D=zT{t!JUE`m_{-}W`GwKrZ&x(E)-i&faG=UvBqx&2yE6D-UzaF) z@|p7WtwzUu$;Tjha6E6{n>Lu|8%8_pmpe@B5JC5!7-XAxI+9iiR=#x6If^!v{L^gP z$e1tr86*#mXW6{>CnorW(SPW;NbP7JLHGG%MX$ho62M+SU=whr8$JDNO?N;q4?uWN0uSeju&b zQXT9k+j0CU;(7ad)5QLOr{CDpOJZJT=UV$yfm)ygY(qLI1)5| zF#da6uJKIPgKIpKJY3_Mo=G0A@l5jk)9sqSOH*bj|9u_p zxO1rPi9{28iR~Y|iPw+9zh8|z@$KU)v0qo6`gXa-_U(TcS4QWNuNzie)tl${7pgUm zyT%qTjlI2Oo*({QvJv|`e(ay`y?zwQix*99=<6s;5bGRNuZO{X9c8#i;y+)v1Z8M8 z^mUZ+dG81Rzn|5d1e7#+T)T8&? z`WsmnQFdKUpOf8NL@B6W&&!=tL`g`Flka-2kV2)a|4*nA<>%}yq=GDaM4hxMqPD$r z>>Zp|L=`WI-Z={3-g`z*+jF)iuITip%vn1MsS$?nPyJz8MD6W-bkw_)B5HS;o^F;` z5v8@^YkS{rHFEVAUn)5sfbdy$z84LQABAr(>f!YT@&}L z<=azI+X^Y2;Nf>!pf9u$n*ey~?5(tv`BoDdpOA zrg3EVQfjco1&^h!HE~;48qF#6FQEq9?exJ;wUoNYJhdjyzij`ZEi+4~3Jq!5*tVtA#P|97 z$Gl3ZfHz$~&XO#p#_J>ukse&9?!__{PL*rU( zxg|5M?Q}rnn&?C{u2~I3)YiT+qX5x!>{kN z;JYXIhb?!Vj^v*ruIhUqo^>Dj1IH8E68>*29;jlbV0A>_~bJS6epYyXp5_WHs~bzgz?6pWUlxNv0Gr>?$V z2fq)x@#t#@czN5}K`CPV z{?6#@uLa)1*F_Yor;y1WPYSO~#cRRnjCguE+1tATjav2y_36cm!E$UcHwPcY8LEqX2$!6D6y zV&fL`l5s2`eWNK}m;O-VbL4`*t3QF4S>CF-9|PxWq$GBWBjXsqdAWO657O}W59X4vuhD)S-E5Asr{>$7YO< z(-<1*ICY5W&geL8U3-whEn;ahe9VeY+q~oN29_cuV z>_9qByG9}%r&pp#$EjQo={SvMb)3eHMmkQ{6OfKm>+wj(DI_VC(QyjDCWUmIh9MoN z!K{wcg-J-qNhJj7I4z1mI!` zCuMD<Nu?|LOM=+vXG9`fmD6gNM#sr?iyx!owC%y~3~nZ?Hqvo2u|PUbZ>*7yQ~SyTjE>XYMFNvTM zKsrwI(vXhR_9Ud^bbLC}ajN=I#ppQc1t44{tK+o$2GVg_Z;f=E+G-*lCkYdz<78Tb zbez&hAe;`XigcV5-B~>*R>$dXG176;$wxX)e08MbwD2_2amp!1I!+-S5iWq$ zaSC-nI!*^Fk&csV2GVgVu|hgdB~y@&(}e*Dm%-{dEggw;oIX89I!>pgkdD*psYu63 zeHPMj+7^m*jw)LqJtkQdq~mn%H>Bfq(Fo}{eVBuEoF<({I!@LyNRMePtK%ep4e2;3 zFGM;{aYK-fQ^;MUgwQZ&+Wnw{`%9j9F4y+M7?<-X0(n*{vc;BVaZ3XA+g9w48~ zmDsk)isLs-Xh-AGaY|(0?ESYoPL72Cf1!~3SDRAT`%3bCZ@v9u9~A#e|7vZ8SIWIL+vQbe!DJUSxEfDho=Hj?>IDwetMQ>NpK~f^?i@b&!tJvg;*`j+4Ws zJ&cZ%*~31Jj?>dQ*|qY_Vs)G*zW;;Kang@OI!oaT!n9j7WAq~kO*7U?*JKAywqIK5kE$>=y8IM%<``V3@soQzH(9j6iYNXJQI zG176mlNiV7IB69YFgi|BchA?Vo3xwNak{C7bev>Ns(okHeP%o>#?I z(D~I8VxE=jd@DKc%60ygoQLH)A4|^5a-E+g=V`gl*OK$LT<34ed0ejXgXeY0I1)5| zFdkd3^OG2l=b5n;biO%(n0Mwn|4hz9bDfVS=cPH%PZRUhoad{Fd26oo*W^4l*ZFL6 zUYl$D;CXH`js%S#jQ`#i$Mc3Dq~jEnfpnbuIv^bz5rW3nDNo*%xt z%jh^gcSbr+iLH^2)1rF{jE>W?W=O{g$1%x+<2mI7(s6QdKsrv|G}3XpzoiAE<5Y0> z6QkpV zJkoJ$9zBfFagwz<$LKiW{0ou?$8+#(q~lbz3+Xtm@<2LH&XbXjlc)mHal&y-^5A&B z8;^9Hq+1~!r-dz%j?-Z|q~r8_(k@2F3CA(XgX4M9Nu=X+T^#8+dA?ZA=s2l1Lpn|) zon;suCo+ynzJ{``#OgSmo`-atLR65B)3TRH$0_>?(s5Gwz+-fr`faaW&nJ`Baf29)y{K=)p4@8j&z(>Ya$(|i^WLCXzN)rn5Ru$8?a6)4(L8 z$ekAf)3YlZ8hm`57lo~$`B7w^6xVzyGH;4&{uG%<#WkOb%&X#> zUq$9wam}|P^RCu$%)cV@u(-w#&dVa>NYMDfcx<`m7hyck6T?=}d@(X_jBEZFnMcMo zpNz~aGal-Aa=8XVG!WCDu9My6!Y{y|EXaByh$*U$+EhX#OlZPVoG5(MO6Ifb@|B zwLP9AeI!9Gj#HsNl29K>sE>rN1|CnSk5ugo;itMi@HQmu(MS~sYKf0uQ7N!f=XAIm z)Z2klS?-yutm;4&$1FaS4Q#7D`djSP9H_Xx$&1C*9lq!#;ins{tHf0w3FFcXC#-i_ z0Q3VV_%)5c97v}+Y?yhkM-Z*$cl3(F+#ouYU#Mqx5o~V|+ZNoF|3&}rgiXufF+g{4 zT*#9`Ctx2O$sVvASkGr$L*syL+G$bgL14p-ry5)b7J7us__nj`e%mo)rs@l759ON7 zIK!K&FQ~<5t}Z`wX*c8~`Iv3?eh04QRP#9v-FbdL9VaDL;RxT+;+BTs4-LxV=YQxM z_?Jat529`IKYm_w9Q5sSB<@$=t{GRlOp|}Sj+3koIx`Jh$0=|qswIx6Up+pA?fcpT zbw2;I>zv}!|CWvuxxOIR=js2~agqX+=K4JS>+3k_vg0D^7xDRR=sHfh!@fPv;r`3l z@Z068+uJxk)f^AMzR$w?P1thRaWd6Io^RN3Rp0yYz>&i7uyI@&wYRda@lc6a$NPJ~ z53lb!{pUJP?mzYOQ0}O2_v>H!f8spVewG{0UcuvvtBzA(9q|>|pq!2MYRGl=tM3=N z>NrU;?>x}x$^$BE=t;hcU7!E&>QX3}3&)=Z6)-=zGI;z+W{>O6+RShi7j@}UByU6C z%d*{MUF!fT??45NcCiInejKmLNoc zJ%vE_c>(vBY{2nC@P`IvNj(KwXA}WA4s_XV3m3Gbro z4wq93MW1N)y0>P_)~eAeT?E*HpdCr(8BaePaRL)$Ej zv#G0JB9B|&bDHL)uj|LTLPFqA=>a@W;&BSceGVJ9xKWKqPr-LGx+4BtJq0C;@6%7< z^Vg`XaQzyZIp0D*2>k#9;RnXS_=S#7`hdKR>HeEGdS3{ldyMKlUb`8eUU7b<_rb}0 zdef9QOUx4Zbe_xkp>HLF>EHqHUabK6Ax7SQ49@YBY4ccsyX^D(@<{-9Z?R`k0l-PO z)=ZoOaOQ(1j8>EiruS4fuY|Z&};HE1W?;Htm$Gldi)r>b5~0B6$q zbUQhKbC8!{aB9O{7~It(uk{@P&ShHQ>s|n-5MU{v1aLDYIyw&nxU4;61D}F=Ud#N* z)bp^_6b9GH?x?ONz`3POnfywePapmBS>!x`+q~$tmOQ}OSk{=zd1u6!=SC8W9J_ygSV?jDE60^D(_ zJfoWcS3aiXt~bDS6CWTW2XNg!=Pa^QYo;?%_NMJ{`FI*FKplC4gDFM z!-meIu7Egs>OP_8Zvc1vXpVavz!{krNWZwvr}xCm^Be*0RM{g2=f78p!O3ODkDUu} zO82$fUm6I?eNR!49@Yvurd>X`!u+JSpvYx zxZkwu1#lM*CKt~IxCaL(>vQzfIoTMGgZS-sEXSF0fRN^VX-p zZjACb74IEL&(5ikR|b|J=qY^z?9b|faAGSKw}!Yl)kX0Tx9DxVy%2ZoV4*d{8Jro{ z6XIrtA>36K=kBsp4&pw!W7Z*hAh4<265** z+=zy_bqA*>K%DwM7jua7*Y43C;=ER&7kxIs8{%51KU@NFZP!F-LR_qwdPOf_oBJYMG!Nm_ z_&v`-oJ32T=@6&zWT!jC?Mm3&1>%}n85BUA#4Lo9Yl?8oo$s7w>M5{4!YK_B)nsr(Mtp?0&nmhdAKVUDIA)LH1!cCo?qYiO5 ztaaZ*+_NO*_7JD)F){|?Hm}m+L7dn)gwq&;aIYV~R)V+=ZO1)_IJNWI(h$e{W6f%a zbBs7u0qdzX3gHxW5w6>TCDIT#Zlz&9#O<+=dH)(%gKG!pL)^}^2L%vkt&ebnR1vQC zkr$QOpq&2-yKIQt`KWaf#Pu1{%K_qU+c!T8adFxR*YY=n<0V)YK-^MM+cb!~W%?l( z;%wfzkAt{BM`-MWxU}8~ccdG_nf1=U2XV&P&r%`I$S{P#rFM~FaOwfMdm+wV9pUPF zE>|Mn8`Sq)p5rKdE@$px8|gptxm=0$gEjl!3A?57=qXgOZ}$FMJ%vod|M)xz?YBKp zSHDQU@2$6AT!Z>WrEecj`0X$10kU7zWZUV>_71q&KT1!*F#^S-P)`AP{s+EiV9vAf zc(l+D!nhM&cOJs)0Cohz{!kbP!Z;9K2S0m10Au(&{VAx-zNi7-+a+ozK1v7*qO)YK z4v5SNq9fLK)au=dPnUM_l^^WPr$0}3GTCsBPy0pfdlmrR_mhVoykQ~3(|&vV2|0kL zpm=kOy?`gKpJefHz>`-zso503qoARieE{%GZh5>`o}4BN7#^+Yd4mrDp3&E4JW2#S zXZr5Q83=gV2J$pa0Z-7;>?uirM>I66R-XPUUJQ@dXU%pS0Z-hrD!r9}=Ur9jeoBC6 z*Zq!TwE@rj4<=%vfJb&(X01F^kJ~dm70>Ks+yT$Ra|#!10Z-5cwFj>|^698{vXWAO zCu&{$O%8xZr+infJV*Q|Fgyo>P6Yf8cr?pg&-MX4XEnEop8z~%XA9=tbmG%T3gc&J z0UonE8*1f==MQIiG_Gvcc^?==zf`|4^i6gU9afQ*5Da(>SIPTt1w6mSE?D*UET8s% zzq(eQ<-Z}GmP1v$6#|~6`%QWj0G{z_PuF+>p2Nddy5=O-~9&og5y=zMdHd1sFE&*VHb*ZF92UYhItG&xVrb-tRMx8^#3 zP0nL;ozEucwYkO*p64dxNYMDf`0s6TJZIYHy-pj*r(5z3qCdOx={Bjw(sPgT>5{2h z`?i97i(R^dTP1i`=TDIJ!1)x%D&1e@1D^H)rF2g>KHYk``T97(vr2j0s55{^HOMXI zBH+PsO!DA(PU4r9W&@teo3kIu0iLlw_Y7A8o~2SE9{T{#>n2Sz_5&Uq$0QGq=WSb; zPYDM+3nvv!EClWP5dAJ^D&SGu_O??1;1N}QT@?s;a2%67IG*h@ck3Gi9=qpCBC&v{ z>fyN63V=t-@C@G&@Jy}jQ*Howa2%67IG)E&9D3ua9?0(*rFnBW;AuN6KJa)J$ghZ9 zwZ6PXFm2jha%-tqbv_2kgX1|+anj-yfX8)SFP9P*KAoQNPHzU_DUOYb2nRd{J+tzH z01u91k_X3g$)86QMggAdNq@Sf0iImt?FM}SkBPL#eJjAzamT*NW`Ku`W0J3-Y|Bh@ zVe%0y)Vmyr2lvr;iay+{%Eh@HF+?5(gJr(XaYxb3E4D=-iKd`N6D+)vj3r5)lK zDLXR<@;r(hQVez1LR56sGY~k81Zawuf*_#$?JOi9pF+8-0Bg3<$Z|yu}9B|#{noom| zhx4Ma6*NDJ%#-4pFGc1}am}A1^QgGyQ;~U9T=T2QJS(pGR%G54*ZeCo4~uL3;Jhp{ zjs%S#jK`L1ei6pwJTYtq%@-r{#<=E>k$GfX^U27(GOqb$WS$wl;ADQcoz360sGyeIyB`8C)q3bEs?4MD@eZO^-%h`SyM>3AI5d5J*Sv>B`vbB*I?4($=*;}L>nLNF`(#*I9#BlB%)5CQoDEH( z<~wW4nevkI~69c{jfTP>i3{U?7^S%0BRP2Ao`CF71Mq*998 zTg`5^GKCt^(ur5GG?|L)<6E-aE1Aj-2)2IGu_o@>q=YS3 zH~Q+dNNQaZ*K&i+V%rv}l<4WPmI@0~D4mEco3GAErlhMjo%M1^rv9+iF&{5e6W4T; z>UmqqOzK4E^%-dv>C~ayT>`GmPp4julE@A3luk`LDi!FxD5J*tOY4!Db+JVz^}%$t zQg47WR*#n(0dRgjd+m|}IID*rt~~UwiQAtue|xLunbdJTaobfB(kWigqKqAL)2UPa zWruimNT|?pLYm-dMqhs%L9Zb`y@O`KDs(GbT`H{~=Tg#_Y2R0exIs4SaW$rDJd;TGV z(z~o@Fmz-(RcSwZipi{WYPz(?0^^qH)Oc$CnDHJpaaHEaJBGi`prXs7HKz?tr#5sN ze|VZxI+diCZ{pEBoqC}gB&+CB6L+_#mrSeY8Po(f-aXsF>C{3yneG<0=~Rqsuao>H z>C{!#$1`$f)WqTINoapA50tH39FnB$KxN44rDOs7`osOpmw`<(o%6B~*e5NjTD<}` zt+zuzF`yG7k#u5QTVTI#WjV23Q)dE}H@&5=Pg<=1L0V>PVQ&KD{_kM;ZUF0!kbHQX z!PkU)-@hF66;QX@Au&KU8nzCM}cPheRTa-xBusO zAF8nD+x_~N{+~Ec)$MGYp9GI9t~v@KZ2wQjS73v3Jl0XL9)sd(5xXvR_5C7O9R;}J z3Xg|dkB^hcu1lYjf{s_;_`}x z5YX1OUXWudVh-H<_u|(gLm0mTy8f3cPR1{O8c1a02RGK;#?8*flgGpNQLZ+go>ulQ zvwYYBPHrxCURH1wbhHDma0Krm4)lv} z)_1>vy@D-%UTiNS!YBtpzoPHZnGoYiFxyRY5<*gJ8V-!AucEgIfdeR1{ApKlWDxc>}kc+LPd z65JnL@5{65@JqwG&VHOfZu}(vsQr0@9G~iYo`^e&`VHKVG{}oVSgC)_19jbx=@R!h z-}`>7zUvu}X4rDSk0@uaXEMSYR{uucM@XE&EwBC5x@> uHXn&=XxsJ0)jz+F;Kd;i;WL0e(Ae)I6zcFx!@AA{j}NZ*RhmDAKmG^db*@AJ diff --git a/test/test_interface.py b/test/test_interface.py index a545106..493be42 100644 --- a/test/test_interface.py +++ b/test/test_interface.py @@ -203,17 +203,23 @@ def test_create_receiver(self): for panel in rec.panels.values(): self.assertEqual(tube_dict_unit["ass_tube_per_panel"], panel.ntubes) - def test_set_tube_pressure_bcs(self): + def test_cycle_tube_pressure_bcs(self): times = np.array([0, 1]) - cyclic_times = np.array([0, 1, 2, 3]) tube_pressure = 12 - pressure = np.ones_like(cyclic_times) + pressure = np.ones_like(times)*tube_pressure + pressure[0] = 0.0 tube = interface.create_tube(tube_dict_unit, times, []) + tube_pressure_bc = receiver.PressureBC(times, pressure) + tube.set_pressure_bc(tube_pressure_bc) + tubes_dict = {"0": tube} - interface.set_tube_pressure_bcs( - tubes_dict, tube_pressure, pressure, cyclic_times + num_cycles = 3 + cyclic_times = np.array([0, 1, 2, 3]) + interface.cycle_tube_pressure_bcs( + tubes_dict, num_cycles, cyclic_times ) - self.assertTrue(np.array_equal(tube_pressure * pressure, tube.pressure_bc.data)) + gold = np.array([0,12,12,12]) + self.assertTrue(np.array_equal(gold, tube.pressure_bc.data)) def test_set_and_downsample_tube_temp_bcs(self): times = np.array([0, 1, 2]) @@ -378,7 +384,7 @@ def test_simple_model(self): # NOTE: making very large for test for single iter convergence pct_err_outlet_temp = 100.25 panel_flow_path = [["1"], ["0"]] - mass_flow_per_path = np.array([650, 650]) # kg/s + mass_flow_per_path = np.ones((len(panel_flow_path),len(times)))*650 # kg/s T_in_per_path = np.array([500, 500]) # Celcius use_cycle_reset_heuristic = False save_heat_to_vtu = False @@ -496,6 +502,8 @@ def test_simple_model(self): rec, tube_dict, manifold_tube_dict, mat_fluid, outlet_p ) inlet_p = outlet_p + interface.convert_Pa_to_MPa(flow_path_p_loss) + interface.update_tube_pressure_bcs(rec, inlet_p, outlet_p) + rec.save(rec_filename+hdf5_ext) # STEP 6: solve structure and life receiver struct_output_dict = { @@ -510,8 +518,6 @@ def test_simple_model(self): analysis_type, is_single_panel_analysis, single_panel_analysis_id, - inlet_p, - outlet_p, num_cycles, solver, struct_output_dict, From 001f74372ff1ca8c0d174684a2faa314749a2366 Mon Sep 17 00:00:00 2001 From: Bryce Mazurowski Date: Thu, 10 Jul 2025 11:33:44 -0500 Subject: [PATCH 6/7] fix: error in SiC.xml data --- srlife/data/damage/SiC.xml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/srlife/data/damage/SiC.xml b/srlife/data/damage/SiC.xml index 64e6f87..f75fc6a 100644 --- a/srlife/data/damage/SiC.xml +++ b/srlife/data/damage/SiC.xml @@ -46,12 +46,12 @@ - 298.15 1073.15 1273.15 1473.15 1673.15 1773.15 - 0.0 0.0 0.0 0.0 0.0 0.0 + 298.15 811.15 1089.15 1366.15 1566.15 + 0.0 0.0 0.0 0.0 0.0 - 298.15 1073.15 1273.15 1473.15 1673.15 1773.15 - 0.0 0.0 0.0 0.0 0.0 0.0 + 298.15 811.15 1089.15 1366.15 1566.15 + 0.0 0.0 0.0 0.0 0.0 298.15 811.15 1089.15 1366.15 1566.15 From 087bcbaa202a90d92d877ec1c191c82832552544 Mon Sep 17 00:00:00 2001 From: Bryce Mazurowski Date: Thu, 10 Jul 2025 13:12:40 -0500 Subject: [PATCH 7/7] fix: formatting fixes to test file --- test/test_interface.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/test/test_interface.py b/test/test_interface.py index 493be42..4b186a6 100644 --- a/test/test_interface.py +++ b/test/test_interface.py @@ -206,7 +206,7 @@ def test_create_receiver(self): def test_cycle_tube_pressure_bcs(self): times = np.array([0, 1]) tube_pressure = 12 - pressure = np.ones_like(times)*tube_pressure + pressure = np.ones_like(times) * tube_pressure pressure[0] = 0.0 tube = interface.create_tube(tube_dict_unit, times, []) tube_pressure_bc = receiver.PressureBC(times, pressure) @@ -215,10 +215,8 @@ def test_cycle_tube_pressure_bcs(self): tubes_dict = {"0": tube} num_cycles = 3 cyclic_times = np.array([0, 1, 2, 3]) - interface.cycle_tube_pressure_bcs( - tubes_dict, num_cycles, cyclic_times - ) - gold = np.array([0,12,12,12]) + interface.cycle_tube_pressure_bcs(tubes_dict, num_cycles, cyclic_times) + gold = np.array([0, 12, 12, 12]) self.assertTrue(np.array_equal(gold, tube.pressure_bc.data)) def test_set_and_downsample_tube_temp_bcs(self): @@ -384,7 +382,7 @@ def test_simple_model(self): # NOTE: making very large for test for single iter convergence pct_err_outlet_temp = 100.25 panel_flow_path = [["1"], ["0"]] - mass_flow_per_path = np.ones((len(panel_flow_path),len(times)))*650 # kg/s + mass_flow_per_path = np.ones((len(panel_flow_path), len(times))) * 650 # kg/s T_in_per_path = np.array([500, 500]) # Celcius use_cycle_reset_heuristic = False save_heat_to_vtu = False @@ -503,7 +501,7 @@ def test_simple_model(self): ) inlet_p = outlet_p + interface.convert_Pa_to_MPa(flow_path_p_loss) interface.update_tube_pressure_bcs(rec, inlet_p, outlet_p) - rec.save(rec_filename+hdf5_ext) + rec.save(rec_filename + hdf5_ext) # STEP 6: solve structure and life receiver struct_output_dict = {