From fe1542cec2a361cbaec0ad8dc019d76c39ac87ac Mon Sep 17 00:00:00 2001 From: Mortara Alessandro Date: Mon, 8 Jan 2024 10:14:38 +0100 Subject: [PATCH 01/26] feat: implemented Class stack, only for some cases The class stack is implemented following the same structure of the strand mix class. The only limit of the new class is related to the evaluations of the self and the mutual inductance, since the old method takes into account the radius of the strand. In this method you must impose a self and mutual inductance. For the future a general way to evaluate these parameters must be implemented. That's the reason why also the conductor.py file has been changed. The new class is able to represent 4 different materials inside the stack. Homogenizes non-superconducting materials into a single stab material and performs the same electrical calculations present in the strand mix class. modified: conductor.py modified: stack_component.py --- source_code/conductor.py | 228 ++++++++-------- source_code/stack_component.py | 458 +++++++++++++++++++-------------- 2 files changed, 380 insertions(+), 306 deletions(-) diff --git a/source_code/conductor.py b/source_code/conductor.py index a25bc5c..0b1e8eb 100644 --- a/source_code/conductor.py +++ b/source_code/conductor.py @@ -3776,91 +3776,92 @@ def __mutual_inductance( """ jj = np.r_[ii + 1 : self.total_elements_current_carriers] - ll = lmod[jj] - mm = lmod[ii] - len_jj = self.total_elements_current_carriers - (ii + 1) - rr = dict( - end_end=np.zeros(len_jj), - end_start=np.zeros(len_jj), - start_start=np.zeros(len_jj), - start_end=np.zeros(len_jj), - ) - - for key in rr.keys(): - rr[key] = self.__vertex_to_vertex_distance(key, ii, jj) - # End for key - - # Additional parameters - alpha2 = ( - rr["start_end"] ** 2 - - rr["start_start"] ** 2 - + rr["end_start"] ** 2 - - rr["end_end"] ** 2 - ) - - cos_eps = np.minimum(np.maximum(alpha2 / (2 * ll * mm), -1.0), 1.0) - sin_eps = np.sin(np.arccos(cos_eps)) + # ll = lmod[jj] + # mm = lmod[ii] + # len_jj = self.total_elements_current_carriers - (ii + 1) + # rr = dict( + # end_end=np.zeros(len_jj), + # end_start=np.zeros(len_jj), + # start_start=np.zeros(len_jj), + # start_end=np.zeros(len_jj), + # ) - dd = 4 * ll ** 2 * mm ** 2 - alpha2 ** 2 - mu = ( - ll - * ( - 2 * mm ** 2 * (rr["end_start"] ** 2 - rr["start_start"] ** 2 - ll ** 2) - + alpha2 * (rr["start_end"] ** 2 - rr["start_start"] ** 2 - mm ** 2) - ) - / dd - ) - nu = ( - mm - * ( - 2 * ll ** 2 * (rr["start_end"] ** 2 - rr["start_start"] ** 2 - mm ** 2) - + alpha2 * (rr["end_start"] ** 2 - rr["start_start"] ** 2 - ll ** 2) - ) - / dd - ) - d2 = rr["start_start"] ** 2 - mu ** 2 - nu ** 2 + 2 * mu * nu * cos_eps + # for key in rr.keys(): + # rr[key] = self.__vertex_to_vertex_distance(key, ii, jj) + # # End for key - # avoid rounding for segments in a plane - d2[d2 < abstol ** 2] = 0 - d0 = np.sqrt(d2) + # # Additional parameters + # alpha2 = ( + # rr["start_end"] ** 2 + # - rr["start_start"] ** 2 + # + rr["end_start"] ** 2 + # - rr["end_end"] ** 2 + # ) - # solid angles - omega = ( - np.arctan( - (d2 * cos_eps + (mu + ll) * (nu + mm) * sin_eps ** 2) - / (d0 * rr["end_end"] * sin_eps) - ) - - np.arctan( - (d2 * cos_eps + (mu + ll) * nu * sin_eps ** 2) - / (d0 * rr["end_start"] * sin_eps) - ) - + np.arctan( - (d2 * cos_eps + mu * nu * sin_eps ** 2) - / (d0 * rr["start_start"] * sin_eps) - ) - - np.arctan( - (d2 * cos_eps + mu * (nu + mm) * sin_eps ** 2) - / (d0 * rr["start_end"] * sin_eps) - ) - ) - omega[d0 == 0.0] = 0.0 + # cos_eps = np.minimum(np.maximum(alpha2 / (2 * ll * mm), -1.0), 1.0) + # sin_eps = np.sin(np.arccos(cos_eps)) + + # dd = 4 * ll ** 2 * mm ** 2 - alpha2 ** 2 + # mu = ( + # ll + # * ( + # 2 * mm ** 2 * (rr["end_start"] ** 2 - rr["start_start"] ** 2 - ll ** 2) + # + alpha2 * (rr["start_end"] ** 2 - rr["start_start"] ** 2 - mm ** 2) + # ) + # / dd + # ) + # nu = ( + # mm + # * ( + # 2 * ll ** 2 * (rr["start_end"] ** 2 - rr["start_start"] ** 2 - mm ** 2) + # + alpha2 * (rr["end_start"] ** 2 - rr["start_start"] ** 2 - ll ** 2) + # ) + # / dd + # ) + # d2 = rr["start_start"] ** 2 - mu ** 2 - nu ** 2 + 2 * mu * nu * cos_eps + + # # avoid rounding for segments in a plane + # d2[d2 < abstol ** 2] = 0 + # d0 = np.sqrt(d2) + + # # solid angles + # omega = ( + # np.arctan( + # (d2 * cos_eps + (mu + ll) * (nu + mm) * sin_eps ** 2) + # / (d0 * rr["end_end"] * sin_eps) + # ) + # - np.arctan( + # (d2 * cos_eps + (mu + ll) * nu * sin_eps ** 2) + # / (d0 * rr["end_start"] * sin_eps) + # ) + # + np.arctan( + # (d2 * cos_eps + mu * nu * sin_eps ** 2) + # / (d0 * rr["start_start"] * sin_eps) + # ) + # - np.arctan( + # (d2 * cos_eps + mu * (nu + mm) * sin_eps ** 2) + # / (d0 * rr["start_end"] * sin_eps) + # ) + # ) + # omega[d0 == 0.0] = 0.0 - # contribution - pp = np.zeros((len_jj, 5), dtype=float) - pp[:, 0] = (ll + mu) * np.arctanh(mm / (rr["end_end"] + rr["end_start"])) - pp[:, 1] = -nu * np.arctanh(ll / (rr["end_start"] + rr["start_start"])) - pp[:, 2] = (mm + nu) * np.arctanh(ll / (rr["end_end"] + rr["start_end"])) - pp[:, 3] = -mu * np.arctanh(mm / (rr["start_start"] + rr["start_end"])) - pp[:, 4] = d0 * omega / sin_eps + # # contribution + # pp = np.zeros((len_jj, 5), dtype=float) + # pp[:, 0] = (ll + mu) * np.arctanh(mm / (rr["end_end"] + rr["end_start"])) + # pp[:, 1] = -nu * np.arctanh(ll / (rr["end_start"] + rr["start_start"])) + # pp[:, 2] = (mm + nu) * np.arctanh(ll / (rr["end_end"] + rr["start_end"])) + # pp[:, 3] = -mu * np.arctanh(mm / (rr["start_start"] + rr["start_end"])) + # pp[:, 4] = d0 * omega / sin_eps - # filter odd cases (e.g. consecutive segments) - pp[np.isnan(pp)] = 0.0 - pp[np.isinf(pp)] = 0.0 + # # filter odd cases (e.g. consecutive segments) + # pp[np.isnan(pp)] = 0.0 + # pp[np.isinf(pp)] = 0.0 # Mutual inductances matrix[ii, jj] = ( - 2 * cos_eps * (pp[:, 0] + pp[:, 1] + pp[:, 2] + pp[:, 3]) - - cos_eps * pp[:, 4] + 5.0e-8 # Mutual inductance imposed by Zappatore input file + # 2 * cos_eps * (pp[:, 0] + pp[:, 1] + pp[:, 2] + pp[:, 3]) + # - cos_eps * pp[:, 4] ) return matrix @@ -3920,23 +3921,24 @@ def __self_inductance_mode1(self, lmod: np.ndarray) -> np.ndarray: for ii, obj in enumerate(self.inventory["StrandComponent"].collection): self_inductance[ii :: self.inventory["StrandComponent"].number] = ( - 2 - * lmod[ii :: self.inventory["StrandComponent"].number] - * ( - np.arcsinh( - lmod[ii :: self.inventory["StrandComponent"].number] - / obj.radius - ) - - np.sqrt( - 1.0 - + ( - obj.radius - / lmod[ii :: self.inventory["StrandComponent"].number] - ) - ** 2 - ) - + obj.radius / lmod[ii :: self.inventory["StrandComponent"].number] - ) + 1.0e-7 + # 2 + # * lmod[ii :: self.inventory["StrandComponent"].number] + # * ( + # np.arcsinh( + # lmod[ii :: self.inventory["StrandComponent"].number] + # / obj.radius + # ) + # - np.sqrt( + # 1.0 + # + ( + # obj.radius + # / lmod[ii :: self.inventory["StrandComponent"].number] + # ) + # ** 2 + # ) + # + obj.radius / lmod[ii :: self.inventory["StrandComponent"].number] + # ) ) return self_inductance @@ -3952,24 +3954,26 @@ def __self_inductance_mode2(self, lmod: np.ndarray) -> np.ndarray: for ii, obj in enumerate(self.inventory["StrandComponent"].collection): - self_inductance[ii :: self.inventory["StrandComponent"].number] = 2 * ( - lmod[ii :: self.inventory["StrandComponent"].number] - * np.log( - ( - lmod[ii :: self.inventory["StrandComponent"].number] - + np.sqrt( - lmod[ii :: self.inventory["StrandComponent"].number] ** 2 - + obj.radius ** 2 - ) - ) - / obj.radius - ) - - np.sqrt( - lmod[ii :: self.inventory["StrandComponent"].number] ** 2 - + obj.radius ** 2 - ) - + lmod[ii :: self.inventory["StrandComponent"].number] / 4 - + obj.radius + self_inductance[ii :: self.inventory["StrandComponent"].number] = ( + 1.0e-7 +# 2 * ( +# lmod[ii :: self.inventory["StrandComponent"].number] +# * np.log( +# ( +# lmod[ii :: self.inventory["StrandComponent"].number] +# + np.sqrt( +# lmod[ii :: self.inventory["StrandComponent"].number] ** 2 +# + obj.radius ** 2 +# ) +# ) +# / obj.radius +# ) +# - np.sqrt( +# lmod[ii :: self.inventory["StrandComponent"].number] ** 2 +# + obj.radius ** 2 +# ) +# + lmod[ii :: self.inventory["StrandComponent"].number] / 4 +# + obj.radius ) return self_inductance diff --git a/source_code/stack_component.py b/source_code/stack_component.py index 4de5212..612042c 100644 --- a/source_code/stack_component.py +++ b/source_code/stack_component.py @@ -78,7 +78,7 @@ cu=density_cu, ge=density_ge, hc276=density_hc276, - re123=density_re123, + ybco=density_re123, sn60pb40=density_sn60pb40, ss=density_ss, ) @@ -89,7 +89,7 @@ cu=thermal_conductivity_cu_nist, ge=thermal_conductivity_ge, hc276=thermal_conductivity_hc276, - re123=thermal_conductivity_re123, + ybco=thermal_conductivity_re123, sn60pb40=thermal_conductivity_sn60pb40, ss=thermal_conductivity_ss, ) @@ -100,7 +100,7 @@ cu=isobaric_specific_heat_cu_nist, ge=isobaric_specific_heat_ge, hc276=isobaric_specific_heat_hc276, - re123=isobaric_specific_heat_re123, + ybco=isobaric_specific_heat_re123, sn60pb40=isobaric_specific_heat_sn60pb40, ss=isobaric_specific_heat_ss, ) @@ -115,6 +115,8 @@ ss=electrical_resistivity_ss, ) +INGEGNERISTIC_MODE = 0 +PHYSICAL_MODE = 1 class StackComponent(StrandComponent): """Class that defines StackComponents objects to model HTS stacks of tapes. @@ -167,6 +169,7 @@ def __init__( current_along=dict(), delta_voltage_along=dict(), linear_power_el_resistance=dict(), + delta_voltage_along_sum=dict(), ) self.dict_scaling_input = dict() # Dictionary initialization: inputs. @@ -190,10 +193,16 @@ def __init__( # Check that costheta is in the range (0,1]. check_costheta(self,dict_file_path["input"],sheet) + self.__compute_cross_section(sheet, dict_file_path["input"]) + self.__get_current_density_cross_section(sheet, dict_file_path["input"]) # Call SolidComponent class constructor to deal with StrandMixedComponent time \ # steps for current, external heating and so on + SolidComponent(simulation, self) + if self.inputs["Stabilizer_material"] != "Cu": + # remove key RRR from inputs if stabilizer is not Cu (cdp, 07/2020) + self.inputs.pop("RRR") if self.operations["IBIFUN"] != -1: # Remove key B_field_units. del self.operations["B_field_units"] @@ -203,39 +212,34 @@ def __init__( # of flag self.operations["IOP_MODE"]. self.deal_with_flag_IOP_MODE() - self.__reorganize_input() - self.__check_consistency(conductor) - # Flag to check if evaluation of homogenized isobaric specific heat can # be done or not (depends on homogenized density evaluation). self.__stack_density_flag = False + self.__reorganize_input() + self.__check_consistency(conductor) # Call method deal_with_fixed_potential to manipulate input about fixed # potential values. self.deal_with_fixed_potential(conductor.inputs["ZLENGTH"]) + + def __repr__(self): + return f"{self.__class__.__name__}(Type: {self.name}, identifier: {self.identifier})" + def __str__(self): + pass + + def __compute_cross_section(self,sheet, file_path:str): + self.cross_section = dict() + self.cross_section["total"] = self.inputs["CROSSECTION"] # Superconductor total cross section in m^2 - self.sc = ( + self.cross_section["sc"] = ( self.inputs["HTS_thickness"] * self.inputs["Stack_width"] - * self.inputs["Tape_number"] + * self.inputs["N_tape"] ) + self.cross_section["sc_sloped"] = self.cross_section["sc"] / self.inputs["COSTETA"] + # Stabilizer (not sc) total cross section in m^2 - self.stabilizer_cross_section = ( - self.tape_thickness_not_sc - * self.inputs["Stack_width"] - * self.inputs["Tape_number"] - ) - - def __repr__(self): - return f"{self.__class__.__name__}(Type: {self.name}, identifier: {self.identifier})" - - def __str__(self): - pass - - def __reorganize_input(self): - """Private method that reorganizes input data stored in dictionary self.inputs to simplify the procedure of properties homogenization.""" - # Create numpy array of string with the identifier of tape material self.tape_material = np.array( [ @@ -258,7 +262,7 @@ def __reorganize_input(self): value.lower() for key, value in self.inputs.items() if key.endswith("material") - and key != "HTS_material" + and key != "superconducting_material" and value != "none" ], dtype=str, @@ -291,6 +295,51 @@ def __reorganize_input(self): ) # Total not superconducting tape thickness in m. self.tape_thickness_not_sc = self.material_thickness_not_sc.sum() + + self.cross_section["stab"] = ( + self.tape_thickness_not_sc + * self.inputs["Stack_width"] + * self.inputs["N_tape"] + ) + + self.cross_section["stab_sloped"] = self.cross_section["stab"] / self.inputs["COSTETA"] + + # Evaluate alpha_0 (stabilizer non stabilizer ratio defined wrt the + # not segregated stabilizer cross section only). + self.cross_section["alpha_0"] = self.cross_section["stab"] / self.cross_section["sc"] + + + def __get_current_density_cross_section(self, sheet, file_path:str): + """Private method that evalutates cross section used to compute the current density and the current sharing temperature according to the definiton of the critical current density scaling parameter c0. + + Args: + sheet (Chartsheet | Worksheet | ReadOnlyWorksheet): sheet with input data for StrandMixedComponent definition. + file_path (str): path to the input file. + + Raises: + ValueError: if self.inputs["C0_MODE"] =! 0 or self.inputs["C0_MODE"] != 1 + """ + + if self.inputs["C0_MODE"] == INGEGNERISTIC_MODE: + # Ingegneristic definition for c0 is used, i.e., the value is + # normalized with respect to the total perpendicular cross section + # of superconducting strand. + + # Convert the ingegneristic definition to the physical definition + # (normalized with respect to the total perpendicular cross section + # of the superconductor) in order to use the superconducting cross + # section in the electric model. + self.inputs["c0"] = self.inputs["c0"] * (1. + self.cross_section["alpha_0"]) + elif (self.inputs["C0_MODE"] != PHYSICAL_MODE + and self.inputs["C0_MODE"] != INGEGNERISTIC_MODE + ): + raise ValueError( + f"Not valid value for flag self.inputs['C0_MODE']. Flag value should be {PHYSICAL_MODE = } or {INGEGNERISTIC_MODE = }, current value is {self.inputs['C0_MODE'] = }.\nPlease check in sheet {sheet} of input file {file_path}.\n" + ) + + + def __reorganize_input(self): + """Private method that reorganizes input data stored in dictionary self.inputs to simplify the procedure of properties homogenization.""" # Create numpy array with density functions according to the tape # material; order is consistent with values in self.tape_material. @@ -436,15 +485,14 @@ def stack_isobaric_specific_heat(self, property: dict) -> np.ndarray: # Set flag to false to trigger error in the next homogenized isobaric # specific heat evaluation if not done properly. self.__stack_density_flag = False - isobaric_specific_heat = np.array( - [ - func(property["temperature"]) - for func in self.isobaric_specific_heat_function - ] + isobaric_specific_heat = np.zeros( + (property["temperature"].size, self.inputs["Material_number"]) ) + for ii, func in enumerate(self.isobaric_specific_heat_function): + isobaric_specific_heat[:, ii] = func(property["temperature"]) # Evaluate homogenized isobaric specific heat of the stack: # cp_eq = sum(s_i*rho_i*cp_i)/sum(s_i*rho_i) - return (isobaric_specific_heat.T * self.__density_numerator).sum( + return (isobaric_specific_heat * self.__density_numerator).sum( axis=1 ) / self.__density_numerator_sum @@ -498,10 +546,8 @@ def stack_electrical_resistivity_not_sc(self, property: dict) -> np.ndarray: electrical_resistivity[:, ii] = func(property["temperature"]) if self.inputs["Material_number"] - 1 > 1: # Evaluate homogenized electrical resistivity of the stack: - # rho_el_eq = s_not_sc * (sum(s_i/rho_el_i))^-1 for any i not sc - return self.tape_thickness_not_sc * np.reciprocal( - (self.material_thickness_not_sc / electrical_resistivity).sum(axis=1) - ) + # rho_el_eq = A_not_sc * (sum(A_i/rho_el_i))^-1 for any i not sc + return self.cross_section["stab_sloped"] * np.reciprocal((self.cross_section["stab_sloped"] / electrical_resistivity).sum(axis=1)) elif self.inputs["Material_number"] - 1 == 1: return electrical_resistivity.reshape(property["temperature"].size) @@ -575,14 +621,22 @@ def solve_current_divider( raise ValueError( f"Arrays rho_el_stabilizer and critical_current must have the same shape.\n {rho_el_stabilizer.shape = };\n{critical_current.shape}.\n" ) - + if rho_el_stabilizer.shape != current.shape: + raise ValueError( + f"Arrays rho_el_stabilizer and current must have the same shape.\n {rho_el_stabilizer.shape = };\n{current.shape}.\n" + ) + if critical_current.shape != current.shape: + raise ValueError( + f"Arrays critical_current and current must have the same shape.\n {critical_current.shape = };\n{current.shape}.\n" + ) + # Evaluate constant value: # psi = rho_el_stab*I_c^n/(E_0*A_stab) psi = ( rho_el_stabilizer * critical_current ** self.inputs["nn"] / self.inputs["E0"] - / self.stabilizer_cross_section + / self.cross_section["stab"] ) # Initialize guess. @@ -714,174 +768,190 @@ def get_electric_resistance(self, conductor: object) -> np.ndarray: np.ndarray: array of electrical resistance in Ohm of length {conductor.grid_input["NELEMS"] = }. """ - critical_current_node = self.sc_cross_section * self.dict_node_pt["J_critical"] critical_current_gauss = ( - self.sc_cross_section * self.dict_Gauss_pt["J_critical"] - ) - - # Get index that correspond to superconducting regime. - ind_sc_node = np.nonzero( - self.dict_node_pt["op_current_sc"] / critical_current_node < 0.95 - )[0] - ind_sc_gauss = np.nonzero( - self.dict_Gauss_pt["op_current_sc"] / critical_current_gauss < 0.95 - )[0] - # Get index that correspond to current sharing regime. - ind_sh_node = np.nonzero( - self.dict_node_pt["op_current_sc"] / critical_current_node >= 0.95 - )[0] - ind_sh_gauss = np.nonzero( - self.dict_Gauss_pt["op_current_sc"] / critical_current_gauss >= 0.95 - )[0] - - # Initialize electric resistance arrays in both nodal and Gauss points; - # this is the equivalent electrical resistance, thus it is defined in - # this way: - # R_eq = R_sc if superconducting regime - # R_eq = R_sc * R_stab/(R_sc + R_stab) is sharing or normal regime - self.dict_node_pt["electric_resistance"] = 10.0 * np.ones( - self.dict_node_pt["temperature"].shape - ) - self.dict_Gauss_pt["electric_resistance"] = 10.0 * np.ones( - self.dict_Gauss_pt["temperature"].shape - ) - - ## SUPERCONDUCTING REGIME ## - - # Strand current in superconducting regime is the one carriend by the - # superconducting material only. - self.dict_node_pt["op_current"][ind_sc_node] = self.dict_node_pt[ - "op_current_sc" - ][ind_sc_node] - self.dict_Gauss_pt["op_current"][ind_sc_gauss] = self.dict_Gauss_pt[ - "op_current_sc" - ][ind_sc_gauss] - - # Initialize array of superconducting electrical resistivit in nodal and Gauss points to None. - self.dict_node_pt["electrical_resistivity_superconductor"] = np.full_like( - self.dict_node_pt["temperature"], None - ) - self.dict_Gauss_pt["electrical_resistivity_superconductor"] = np.full_like( - self.dict_Gauss_pt["temperature"], None - ) - - # Compute superconducting electrical resistivity only in index for - # which the superconducting regime is guaranteed, using the power low. - self.dict_node_pt["electrical_resistivity_superconductor"][ - ind_sc_node - ] = self.superconductor_power_law( - self.dict_node_pt["op_current_sc"][ind_sc_node], - critical_current_node[ind_sc_node], - self.dict_node_pt["J_critical"][ind_sc_node], - ) - self.dict_Gauss_pt["electrical_resistivity_superconductor"][ - ind_sc_gauss - ] = self.superconductor_power_law( - self.dict_Gauss_pt["op_current_sc"][ind_sc_gauss], - critical_current_gauss[ind_sc_gauss], - self.dict_Gauss_pt["J_critical"][ind_sc_gauss], - ) - # Evaluate electic resistance in superconducting region (superconductor - # only). - self.dict_Gauss_pt["electric_resistance"][ - ind_sc_gauss - ] = self.electric_resistance( - conductor, "electrical_resistivity_superconductor", ind_sc_gauss - ) - - ## SHARING OR NORMAL REGIME ## - - # Strand current in sharing regime is the one carried by the both the - # superconducting and the stabilizer materials. - # self.dict_node_pt["op_current"][ind_sh_node] = self.dict_node_pt["op_current"][ind_sh_node] - # self.dict_Gauss_pt["op_current"][ind_sh_node] = self.dict_node_pt["op_current"][ind_sh_gauss] - - # Evaluate how the current is distributed solving the current divider - # problem in both nodal and Gauss points. - sc_current_node, stab_current_node = self.solve_current_divider( - self.dict_node_pt["electrical_resistivity_stabilizer"][ind_sh_node], - critical_current_node[ind_sh_node], - self.dict_node_pt["op_current"][ind_sh_node], - ) - sc_current_gauss, stab_current_gauss = self.solve_current_divider( - self.dict_Gauss_pt["electrical_resistivity_stabilizer"][ind_sh_gauss], - critical_current_gauss[ind_sh_gauss], - self.dict_Gauss_pt["op_current"][ind_sh_gauss], + self.cross_section["sc"] * self.dict_Gauss_pt["J_critical"] ) - # Get index of the normal region, to avoid division by 0 in evaluation - # of sc electrical resistivity with the power law. - ind_normal_node = np.nonzero( - ( - stab_current_node / self.dict_node_pt["op_current"][ind_sh_node] - > 0.999999 + # Make initialization only once for each conductor object. + if conductor.cond_num_step == 0: + # Initialize electric resistance arrays in Gauss point; this is the + # equivalent electrical resistance, thus it is defined in this way: + # R_eq = R_sc if superconducting regime + # R_eq = R_sc * R_stab/(R_sc + R_stab) is normal regime. + self.dict_Gauss_pt["electric_resistance"] = np.full_like( + self.dict_Gauss_pt["temperature"], None ) - | (sc_current_node < 1.0) - )[0] - ind_normal_gauss = np.nonzero( - ( - stab_current_gauss / self.dict_Gauss_pt["op_current"][ind_sh_gauss] - > 0.999999 + + # Initialize array of superconducting electrical resistivit in + # Gauss point to None. + self.dict_Gauss_pt["electrical_resistivity_superconductor"] = np.full_like( + self.dict_Gauss_pt["temperature"], None ) - | (sc_current_gauss < 1.0) - )[0] - - ## NORMAL REGIME ONLY ## - if ind_normal_node.any(): - # Get the index of location of true current sharing region; - # overwrite ind_sh_node. - ind_sh_node = np.nonzero( - ( - stab_current_node / self.dict_node_pt["op_current"][ind_sh_node] - <= 0.999999 - ) - | (sc_current_node >= 1.0) - )[0] - if ind_normal_gauss.any(): - # Get the index of location of true current sharing region; - # overwrite ind_sh_gauss. - ind_sh_gauss = np.nonzero( - ( - stab_current_gauss / self.dict_Gauss_pt["op_current"][ind_sh_gauss] - <= 0.999999 - ) - | (sc_current_gauss >= 1.0) - )[0] + + # Get index for which abs(critical_current_gauss) == 0 (inside normal + # zone by definition). + ind_zero = np.nonzero(abs(critical_current_gauss) == 0)[0] + # Get index for which abs(critical_current_gauss) > 0 (outside normal + # zone by definition). + ind_not_zero = np.nonzero(abs(critical_current_gauss) > 0)[0] + + # Check if np array ind_zero is not empty: NORMAL REGION BY DEFINITION + if ind_zero.any(): # Evaluate electic resistance in normal region (stabilizer only). self.dict_Gauss_pt["electric_resistance"][ - ind_normal_gauss + ind_zero ] = self.electric_resistance( - conductor, "electrical_resistivity_stabilizer", ind_normal_gauss + conductor, "electrical_resistivity_stabilizer", "stab", ind_zero ) - ## SHARING REGIME ONLY ## - # Evaluate the electrical resistivity of the superconductor according - # to the power low in both nodal and Gauss points in Ohm*m. - self.dict_node_pt["electrical_resistivity_superconductor"][ - ind_sh_node - ] = self.superconductor_power_law( - sc_current_node[ind_sh_node], - critical_current_node[ind_sh_node], - self.dict_node_pt["J_critical"][ind_sh_node], - ) - self.dict_Gauss_pt["electrical_resistivity_superconductor"][ - ind_sh_gauss - ] = self.superconductor_power_law( - sc_current_gauss[ind_sh_gauss], - critical_current_gauss[ind_sh_gauss], - self.dict_Gauss_pt["J_critical"][ind_sh_gauss], - ) + # Check if np array ind_not_zero is not empty: deal with index that + # are outside normal zone by definition; however some of them could + # still identify a normal region. + if ind_not_zero.any(): - # Evaluate the equivalent electric resistance in Ohm. - self.dict_Gauss_pt["electric_resistance"][ - ind_sh_gauss - ] = self.parallel_electric_resistance( - conductor, - [ - "electrical_resistivity_superconductor", - "electrical_resistivity_stabilizer", - ], - ind_sh_gauss, - ) + # Get index that correspond to superconducting regime. + ind_sc_gauss = np.nonzero( + self.dict_Gauss_pt["op_current_sc"][ind_not_zero] / critical_current_gauss[ind_not_zero] < 0.95 + )[0] + # Get index that correspond to the normal regime. + ind_normal_gauss = np.nonzero( + self.dict_Gauss_pt["op_current_sc"][ind_not_zero] / critical_current_gauss[ind_not_zero] >= 0.95 + )[0] + + ## SUPERCONDUCTING REGIME ## + + # Check if np array ind_sc_gauss is not empty. + if ind_sc_gauss.any(): + # Current in superconducting regime is the carried by + # superconducting material only. + self.dict_Gauss_pt["op_current"][ind_not_zero[ind_sc_gauss]] = self.dict_Gauss_pt[ + "op_current_sc" + ][ind_not_zero[ind_sc_gauss]] + + # Compute superconducting electrical resistivity only in index + # for which the superconducting regime is guaranteed, using the + # power low. + self.dict_Gauss_pt["electrical_resistivity_superconductor"][ind_not_zero[ind_sc_gauss]] = self.superconductor_power_law( + self.dict_Gauss_pt["op_current"][ind_not_zero[ind_sc_gauss]], + critical_current_gauss[ind_not_zero[ind_sc_gauss]], + self.dict_Gauss_pt["J_critical"][ind_not_zero[ind_sc_gauss]] + ) + + # Evaluate electic resistance in superconducting region + # (superconductor only). + self.dict_Gauss_pt["electric_resistance"][ind_not_zero[ + ind_sc_gauss + ]] = self.electric_resistance( + conductor, "electrical_resistivity_superconductor", "sc", ind_not_zero[ind_sc_gauss] + ) + + ## NORMAL REGIME ## + + # Check if np array ind_normal_gauss is not empty. + if ind_normal_gauss.any(): + + # Evaluate how the current is distributed solving the current + # divider problem in Gauss point. + sc_current_gauss, stab_current_gauss = self.solve_current_divider( + self.dict_Gauss_pt["electrical_resistivity_stabilizer"][ind_not_zero[ind_normal_gauss]], + critical_current_gauss[ind_not_zero[ind_normal_gauss]], + self.dict_Gauss_pt["op_current"][ind_not_zero[ind_normal_gauss]] + ) + + # Get index of the normal region where all the current is + # carried by the stabilizer, to avoid division by 0 in + # evaluation of superconducting electrical resistivity with the + # power law. + ind_stab_gauss = np.nonzero( + ( + stab_current_gauss / self.dict_Gauss_pt["op_current"][ind_not_zero[ind_normal_gauss]] + > 0.999999 + ) + | (sc_current_gauss < 1.0) + )[0] + + ## CURRENT CARRIED BY THE STABILIZER ## + # Check if np array ind_stab_gauss is not empty. + if ind_stab_gauss.any(): + # Get the index of location of current sharing region, if + # any. + ind_sh_gauss = np.nonzero( + ( + stab_current_gauss + / self.dict_Gauss_pt["op_current"][ind_not_zero[ind_normal_gauss]] + <= 0.999999 + ) + | (sc_current_gauss >= 1.0) + )[0] + + # Check if nparray ind_sh_gauss is not empty. + if ind_sh_gauss.any(): + # ind_sh_gauss is not empty. + # Get final index of the location of the current + # sharing zone, keeping into account that + # sc_current_gauss is already filtered on + # ind_not_zero[ind_normal_gauss]. + ind_shf_gauss = {1:ind_sh_gauss,2:ind_not_zero[ind_normal_gauss[ind_sh_gauss]]} + else: + # ind_sh_gauss is empty. + # Get final index of the location of the current + # sharing zone, keeping into account that + # sc_current_gauss is already filtered on + # ind_not_zero[ind_normal_gauss] and that ind_sh_gauss + # is empty. + ind_shf_gauss = {1:ind_sh_gauss,2:ind_not_zero[ind_normal_gauss]} + + # Evaluate electic resistance in normal region (stabilizer + # only). + self.dict_Gauss_pt["electric_resistance"][ + ind_shf_gauss[2] + ] = self.electric_resistance( + conductor, "electrical_resistivity_stabilizer", "stab", ind_shf_gauss[2] + ) + else: + # Get final index of the location of the current sharing + # zone, keeping into account that sc_current_gauss is + # already filtered on ind_not_zero[ind_normal_gauss]. In + # this case the current is shared between the + # superconductor and the stabilizer, so we need to exploit + # all the index in array ind_normal_gauss for the array + # sc_current_gauss. + ind_shf_gauss = {1:np.nonzero(ind_normal_gauss>=0)[0],2:ind_not_zero[ind_normal_gauss]} + + ## CURRENT SHARED BY THE SUPERCONDUCTOR AND THE STABILIZER ## + if ind_shf_gauss[1].any(): + # Evaluate the electrical resistivity of the superconductor + # according to the power low in Gauss point in Ohm*m. + self.dict_Gauss_pt["electrical_resistivity_superconductor"][ + ind_shf_gauss[2] + ] = self.superconductor_power_law( + sc_current_gauss[ind_shf_gauss[1]], + critical_current_gauss[ind_shf_gauss[2]], + self.dict_Gauss_pt["J_critical"][ind_shf_gauss[2]], + ) + + # Evaluate the equivalent electric resistance in Ohm. + self.dict_Gauss_pt["electric_resistance"][ + ind_shf_gauss[2] + ] = self.parallel_electric_resistance( + conductor, + [ + "electrical_resistivity_superconductor", + "electrical_resistivity_stabilizer", + ],["sc","stab"], + ind_shf_gauss[2], + ) + + # Compute voltage along stabilizer. + v_stab = self.dict_Gauss_pt["electrical_resistivity_stabilizer"][ + ind_shf_gauss[2] + ] * stab_current_gauss[ind_shf_gauss[1]] / self.cross_section["stab"] + # Compute voltage along superconductor + v_sc = self.inputs["E0"] * (sc_current_gauss[ind_shf_gauss[1]] / critical_current_gauss[ind_shf_gauss[2]]) ** self.inputs["nn"] + # Check that the voltage along stabilizer is equal to the + # voltage along superconductor (i.e, check the reliability + # of the current divider). + if all(np.isclose(v_stab,v_sc)) == False: + raise ValueError(f"Voltage difference along superconductor and stabilizer must be the same.") return self.dict_Gauss_pt["electric_resistance"] From 373f5329556a31a1c36cee4008e6f12b59d67925 Mon Sep 17 00:00:00 2001 From: Mortara Alessandro Date: Tue, 30 Jan 2024 16:54:38 +0100 Subject: [PATCH 02/26] bug fix: modified the lower limit for T_cs In function current_sharing_temperature_re123 bisection method is used to evaluate the current sharing temperature of the superconductor. In specific cases occurs the error: "r = _zeros._bisect(f, a, b, xtol, rtol, maxiter, args, full_output, disp) ValueError: f(a) and f(b) must have different sign" For this reson T_lower is changed from 3.5 to 0.5, to avoid this problem. modified: properties_of_materials/rare_earth_123.py --- source_code/properties_of_materials/rare_earth_123.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source_code/properties_of_materials/rare_earth_123.py b/source_code/properties_of_materials/rare_earth_123.py index 9722c3b..48aa06f 100644 --- a/source_code/properties_of_materials/rare_earth_123.py +++ b/source_code/properties_of_materials/rare_earth_123.py @@ -393,7 +393,7 @@ def critical_current_density_bisection_re123(TT, BB, JOP, TC0M, BC20M, C): for _, vv in enumerate(JC_ind): - T_lower = np.array([3.5]) + T_lower = np.array([0.5]) # T_upper = np.array([40.0]) T_upper = np.array([TC0M]) From 5a75ebf4baed14cb89c8412ab5425f9791945334 Mon Sep 17 00:00:00 2001 From: Mortara Alessandro Date: Mon, 5 Feb 2024 11:04:56 +0100 Subject: [PATCH 03/26] fix: error in T_cs of re_123 Inside properties_of_materials/rare_earth_123.py the formula: TCSRE123[vv] = optimize.bisect( critical_current_density_bisection_re123, T_lower, T_upper, ex_args, xtol=1e-5, ) The formula gave errors related to the bisection formula. In the interval [a, b] as it was chosen f(a) and f(b) had the same sign. For this reason T_lower, i.e. a, was lowered to 0, so as to solve the problem on the bisection formula. modified: properties_of_materials/rare_earth_123.py --- source_code/properties_of_materials/rare_earth_123.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source_code/properties_of_materials/rare_earth_123.py b/source_code/properties_of_materials/rare_earth_123.py index 48aa06f..21ccf45 100644 --- a/source_code/properties_of_materials/rare_earth_123.py +++ b/source_code/properties_of_materials/rare_earth_123.py @@ -393,7 +393,7 @@ def critical_current_density_bisection_re123(TT, BB, JOP, TC0M, BC20M, C): for _, vv in enumerate(JC_ind): - T_lower = np.array([0.5]) + T_lower = np.array([0]) # T_upper = np.array([40.0]) T_upper = np.array([TC0M]) From 61e4a81dc34dd7bc30db5aea0c5498cf67b48a62 Mon Sep 17 00:00:00 2001 From: DanielePlacido Date: Fri, 9 Feb 2024 11:08:34 +0100 Subject: [PATCH 04/26] chore: remove mod1 to compute Joule power along MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit During the fix of the evaluation of the Joule power along the current carrier, two different approaches were used: * mod1 -> P_joule = R * I^2 * mod2 -> P_joule = Delta_Phi * I The approaches should be equivalent if the inducance is negligible, while mod2 is more general and conservative if the inductance is no longer negligible. After discussion with prof. Zach Hartwig and Dr. Nicolò Riva, we agreed to remove mod1 and keep mod2 for this evaluation. List of changes: * removed key "integral_power_el_res_mod1" from dictionary obj.dict_Gauss_pt of current carrier objects. * renamed key "integral_power_el_res_mod2" into "integral_power_el_res" in dictionary obj.dict_Gauss_pt of current carrier objects. * removed equivalence checks * updated related comments * updated docstring LList of modified methods in class SolidComponents: * initialize_electric_quantities * set_power_array_to_zeros * get_joule_power_along * get_joule_power_along_steady List of modified methods in class Conductor: * __electric_method_steady * __build_heat_source_gauss_pt * eval_integral_joule_power Class: SolidComponents, Conductor modified: conductor.py modified: solid_component.py --- source_code/conductor.py | 38 ++++++++++---------------------- source_code/solid_component.py | 40 +++++++++++++--------------------- 2 files changed, 27 insertions(+), 51 deletions(-) diff --git a/source_code/conductor.py b/source_code/conductor.py index 90973c3..941952e 100644 --- a/source_code/conductor.py +++ b/source_code/conductor.py @@ -4368,8 +4368,7 @@ def __electric_method_steady(self): # conductor_initialization); # 2) all the Joule power contribution due to the electric solution # are evaluated calling method self.__get_heat_source_em_steady. - # N.B. arrays obj.dict_Gauss_pt["integral_power_el_res_mod1"] - # obj.dict_Gauss_pt["integral_power_el_res_mod2"] and + # N.B. arrays obj.dict_Gauss_pt["integral_power_el_res"] and # obj.dict_node_pt["integral_power_el_cond"] are set to zero for # the next evaluation inside method __build_heat_source_gauss_pt self.__build_heat_source_gauss_pt() @@ -4627,8 +4626,7 @@ def __build_heat_source_gauss_pt(self): """Private method that builds heat source therms in Gauss points for strand and jacket objects. Sets to zeros the following arrays as preliminary step for the next evaluation: - * strand.dict_Gauss_pt["integral_power_el_res_mod1"] - * strand.dict_Gauss_pt["integral_power_el_res_mod2"] + * strand.dict_Gauss_pt["integral_power_el_res"] * strand.dict_node_pt["integral_power_el_cond"] """ @@ -4649,8 +4647,7 @@ def __build_heat_source_gauss_pt(self): + strand.dict_Gauss_pt["linear_power_el_resistance"] ) - # Set arrays strand.dict_Gauss_pt["integral_power_el_res_mod1"], - # strand.dict_Gauss_pt["integral_power_el_res_mod2"] and + # Set arrays strand.dict_Gauss_pt["integral_power_el_res"] and # strand.dict_node_pt["integral_power_el_cond"] to zero for the # next evaluation. strand.set_power_array_to_zeros(self) @@ -4700,13 +4697,12 @@ def __build_heat_source_gauss_pt(self): def eval_integral_joule_power(self): - """Method that evaluates the numerator of the expression used to evaluate the integral value of the Joule power. For the Joule power due to the electric resistance along the current carriers two different but equivalent approaches are used: - * mode1 -> P_Joule = R * I^2 (R electric resistance in Ohm, I electric current) - * mode2 -> P_Joule = \Delta_Phi * I + """Method that evaluates the numerator of the expression used to evaluate the integral value of the Joule power. The Joule power along the current carriers is computed as: + P_Joule = \Delta_Phi * I with - * R electric resistance in Ohm * I electric current in A * \Delta_Ph electric voltage potential difference along current carriers in V + This approach was discussed with prof. Zach Hartwig and Dr. Nicolò Riva and is more general and coservative that the evaluation that takes into account only the contribution of the electric resistance. For the Joule power due to the electric conductances between current carriers it is exploited the power computed in method get_total_joule_power_electric_conductance Regardless of the kind of Joule power, the integration can be performed as follows: P_Joule = 1/Delta_t_TH * int_0^Delta_t_TH (dt P_{Joule,i}) @@ -4730,27 +4726,17 @@ def eval_integral_joule_power(self): # Loop on StrandComponent objects. for strand in self.inventory["StrandComponent"].collection: - # Compute the numerator of the integral Joule power due to the - # electric resistance along the current carrier, mode1. - strand.dict_Gauss_pt["integral_power_el_res_mod1"] += ( - strand.dict_Gauss_pt["current_along"] ** 2 - * strand.dict_Gauss_pt["electric_resistance"] - * self.electric_time_step - ) - - # Compute the numerator of the integral Joule power due to the - # electric resistance along the current carrier, mode2. - strand.dict_Gauss_pt["integral_power_el_res_mod2"] += ( + # Compute the numerator of the integral Joule power along the + # current carrier. + # N.B. this evaluation accounts aslo for the voltage due to the + # inductance and is a conservative an more general approach. + # Discussed with prof. Zach Hartwig and Dr. Nicolò Riva. + strand.dict_Gauss_pt["integral_power_el_res"] += ( strand.dict_Gauss_pt["current_along"] * strand.dict_Gauss_pt["delta_voltage_along"] * self.electric_time_step ) - # Check equivalence of Mode 1 and Mode 2 (they should be equivalent - # but Mode 2 may give numerical cancellation). - if not np.allclose(strand.dict_Gauss_pt["integral_power_el_res_mod1"],strand.dict_Gauss_pt["integral_power_el_res_mod2"]): - warnings.warn("P_Joule = R I^2 dt_em != Delta_Phi I dt_em.Possible violation of the energy conservation!") - # Compute the numerator of the integral Joule power due to the # electric conductance between current carriers. strand.dict_node_pt["integral_power_el_cond"] += ( diff --git a/source_code/solid_component.py b/source_code/solid_component.py index e3d98f1..07bbdb2 100644 --- a/source_code/solid_component.py +++ b/source_code/solid_component.py @@ -659,8 +659,7 @@ def initialize_electric_quantities(self, conductor): * self.dict_Gauss_pt["current_along"] * self.dict_Gauss_pt["delta_voltage_along"] * self.dict_Gauss_pt["delta_voltage_along_sum"] - * self.dict_Gauss_pt["integral_power_el_res_mod1"] - * self.dict_Gauss_pt["integral_power_el_res_mod2"] + * self.dict_Gauss_pt["integral_power_el_res"] * self.dict_Gauss_pt["linear_power_el_resistance"] * self.dict_node_pt["total_linear_power_el_cond"] * self.dict_node_pt["total_power_el_cond"] @@ -675,12 +674,10 @@ def initialize_electric_quantities(self, conductor): self.dict_Gauss_pt["current_along"] = np.zeros(n_elem) self.dict_Gauss_pt["delta_voltage_along"] = np.zeros(n_elem) self.dict_Gauss_pt["delta_voltage_along_sum"] = np.zeros(n_elem) - # Introduced to evaluate the integral value Joule power due to electric - # resistance along the cable in the electric module as R * I^2. - self.dict_Gauss_pt["integral_power_el_res_mod1"] = np.zeros(n_elem) + # Introduced to evaluate the integral value Joule power due to electric # resistance along the cable in the electric module as Delta V * I. - self.dict_Gauss_pt["integral_power_el_res_mod2"] = np.zeros(n_elem) + self.dict_Gauss_pt["integral_power_el_res"] = np.zeros(n_elem) self.dict_node_pt["total_power_el_cond"] = np.zeros(n_nod) # Introduced to evaluate the integral value Joule power due to electric @@ -716,12 +713,11 @@ def set_power_array_to_zeros(self, conductor): n_elem = conductor.grid_input["NELEMS"] # Set the following arrays to zeros - self.dict_Gauss_pt["integral_power_el_res_mod1"] = np.zeros(n_elem) - self.dict_Gauss_pt["integral_power_el_res_mod2"] = np.zeros(n_elem) + self.dict_Gauss_pt["integral_power_el_res"] = np.zeros(n_elem) self.dict_node_pt["integral_power_el_cond"] = np.zeros(n_nod) def get_joule_power_along(self, conductor: object): - """Method that evaluate the contribution to the integral power in the element of Joule power (in W/m) due to the electic resistances along the SolidComponent objects. + """Method that evaluate the contribution to the integral power in the element of Joule power (in W/m) due to the voltage difference (due to electric resistance and inductance) along the SolidComponent objects. This approach was discussed with prof. Zach Hartwig and Dr. Nicolò Riva and is more general and coservative that the evaluation that takes into account only the contribution of the electric resistance. This method should be called in the electric method, when the transient solution is used and only for current carriers. Args: @@ -733,14 +729,14 @@ def get_joule_power_along(self, conductor: object): # Finalize the evaluation of the integral value of the Joule power due # to electric resistance along the current carrier started in method - # conductor.eval_integral_joule_power. Array integral_power_el_res_mod1 + # conductor.eval_integral_joule_power. Array integral_power_el_res # stores the numerator (energy in J) that here is divided by the # thermal hydraulic to get again a power (W), which is further divided # by the length of the discreticazion element corrected by cos(theta) # to get a linear power density (W/m): # P_Joule = sum_1^N_em P_{Joule,i} / (Delta_t_TH * Delta_z * costheta) integral_j_pow_along = ( - self.dict_Gauss_pt["integral_power_el_res_mod1"] + self.dict_Gauss_pt["integral_power_el_res"] / (conductor.grid_features["delta_z"] * self.inputs["COSTETA"] * conductor.time_step) ) @@ -823,7 +819,8 @@ def get_joule_power_across(self, conductor: object): ) def get_joule_power_along_steady(self, conductor: object): - """Method that evaluate the contribution to the total power in the element of Joule power (in W/m) due to the electic resistances along the SolidComponent objects. + """Method that evaluate the contribution to the total power in the element of Joule power (in W/m) due to the voltage difference (due to electric resistance and inductance) along the SolidComponent objects. + This approach was discussed with prof. Zach Hartwig and Dr. Nicolò Riva and is more general and coservative that the evaluation that takes into account only the contribution of the electric resistance. This method should be called in the electric method, when the transient solution is used and for the solid component that actually carries a current. It works when the steady state solution for the electric module is computed. Args: @@ -831,31 +828,24 @@ def get_joule_power_along_steady(self, conductor: object): """ # Alias - el_res = self.dict_Gauss_pt["electric_resistance"] current = self.dict_Gauss_pt["current_along"] voltage = self.dict_Gauss_pt["delta_voltage_along"] d_z_tilde = conductor.grid_features["delta_z"] * self.inputs["COSTETA"] - # Mode 1: evaluate Joule linear power along the strand in W, due - # to electric resistances only for current carriers: - # P_along = R_along * I_along ^2 - self.dict_Gauss_pt["integral_power_el_res_mod1"] = current ** 2 * el_res - - # Mode 2: evaluate Joule linear power along the strand in W, due + # Evaluate Joule linear power along the strand in W, due # to electric resistances only for current carriers: # P_along = Delta_Phi_along * I_along - self.dict_Gauss_pt["integral_power_el_res_mod2"] = voltage * current + # N.B. this evaluation accounts aslo for the voltage due to the + # inductance and is a conservative an more general approach. Discussed + # with prof. Zach Hartwig and Dr. Nicolò Riva. + self.dict_Gauss_pt["integral_power_el_res"] = voltage * current - # Check equivalence of Mode 1 and Mode 2 (they should be equivalent - # but Mode 2 may give numerical cancellation). - if not np.allclose(self.dict_Gauss_pt["integral_power_el_res_mod1"],self.dict_Gauss_pt["integral_power_el_res_mod2"]): - warnings.warn("P_Joule = R I^2 != Delta_Phi I. Possible violation of the energy conservation!") # Convert W in W/m keping into account the cos(theta). # This is independent of the time integration method since at the # initialization (when this function is used) BE, CN and AM4 should # fill only the 0 index column. self.dict_Gauss_pt["linear_power_el_resistance"][:, 0] = ( - self.dict_Gauss_pt["integral_power_el_res_mod1"] / d_z_tilde + self.dict_Gauss_pt["integral_power_el_res"] / d_z_tilde ) def get_joule_power_across_steady(self, conductor: object): From 3f597445b592c4d35497ed8f99bd626e35be728f Mon Sep 17 00:00:00 2001 From: Mortara Alessandro Date: Mon, 12 Feb 2024 10:41:48 +0100 Subject: [PATCH 05/26] refactor(indices): fun get_electric_resistance Improved the readability of the code and increased its simplicity by removing the multiple indices and distinguishing only two cases: that in which the critical current is =0 and that in which the critical current is other than 0. In the first case, the equivalent resistance will be set equal to the equivalent resistance of the stabilizer part. In the second case, the current divider will be calculated to derive the current on the superconductor, from it with the power law the resistivity of the superconductor will be estimated, and then the equivalent resistance of the stabilizer and superconductor together will be calculated function: get_electric_resistance modified: strand_mixed_component.py --- source_code/stack_component.py | 181 +++++++-------------------------- 1 file changed, 36 insertions(+), 145 deletions(-) diff --git a/source_code/stack_component.py b/source_code/stack_component.py index 612042c..2fd72e8 100644 --- a/source_code/stack_component.py +++ b/source_code/stack_component.py @@ -772,7 +772,7 @@ def get_electric_resistance(self, conductor: object) -> np.ndarray: self.cross_section["sc"] * self.dict_Gauss_pt["J_critical"] ) - # Make initialization only once for each conductor object. + # Make initialization only once for each conductor object. if conductor.cond_num_step == 0: # Initialize electric resistance arrays in Gauss point; this is the # equivalent electrical resistance, thus it is defined in this way: @@ -787,7 +787,7 @@ def get_electric_resistance(self, conductor: object) -> np.ndarray: self.dict_Gauss_pt["electrical_resistivity_superconductor"] = np.full_like( self.dict_Gauss_pt["temperature"], None ) - + # Get index for which abs(critical_current_gauss) == 0 (inside normal # zone by definition). ind_zero = np.nonzero(abs(critical_current_gauss) == 0)[0] @@ -809,149 +809,40 @@ def get_electric_resistance(self, conductor: object) -> np.ndarray: # still identify a normal region. if ind_not_zero.any(): - # Get index that correspond to superconducting regime. - ind_sc_gauss = np.nonzero( - self.dict_Gauss_pt["op_current_sc"][ind_not_zero] / critical_current_gauss[ind_not_zero] < 0.95 - )[0] - # Get index that correspond to the normal regime. - ind_normal_gauss = np.nonzero( - self.dict_Gauss_pt["op_current_sc"][ind_not_zero] / critical_current_gauss[ind_not_zero] >= 0.95 - )[0] - - ## SUPERCONDUCTING REGIME ## - - # Check if np array ind_sc_gauss is not empty. - if ind_sc_gauss.any(): - # Current in superconducting regime is the carried by - # superconducting material only. - self.dict_Gauss_pt["op_current"][ind_not_zero[ind_sc_gauss]] = self.dict_Gauss_pt[ - "op_current_sc" - ][ind_not_zero[ind_sc_gauss]] - - # Compute superconducting electrical resistivity only in index - # for which the superconducting regime is guaranteed, using the - # power low. - self.dict_Gauss_pt["electrical_resistivity_superconductor"][ind_not_zero[ind_sc_gauss]] = self.superconductor_power_law( - self.dict_Gauss_pt["op_current"][ind_not_zero[ind_sc_gauss]], - critical_current_gauss[ind_not_zero[ind_sc_gauss]], - self.dict_Gauss_pt["J_critical"][ind_not_zero[ind_sc_gauss]] + sc_current_gauss, stab_current_gauss = self.solve_current_divider( + self.dict_Gauss_pt["electrical_resistivity_stabilizer"][ind_not_zero], + critical_current_gauss[ind_not_zero], + self.dict_Gauss_pt["op_current"][ind_not_zero] ) - - # Evaluate electic resistance in superconducting region - # (superconductor only). - self.dict_Gauss_pt["electric_resistance"][ind_not_zero[ - ind_sc_gauss - ]] = self.electric_resistance( - conductor, "electrical_resistivity_superconductor", "sc", ind_not_zero[ind_sc_gauss] - ) - - ## NORMAL REGIME ## - - # Check if np array ind_normal_gauss is not empty. - if ind_normal_gauss.any(): - - # Evaluate how the current is distributed solving the current - # divider problem in Gauss point. - sc_current_gauss, stab_current_gauss = self.solve_current_divider( - self.dict_Gauss_pt["electrical_resistivity_stabilizer"][ind_not_zero[ind_normal_gauss]], - critical_current_gauss[ind_not_zero[ind_normal_gauss]], - self.dict_Gauss_pt["op_current"][ind_not_zero[ind_normal_gauss]] + + self.dict_Gauss_pt["electrical_resistivity_superconductor"][ + ind_not_zero + ] = self.superconductor_power_law( + sc_current_gauss, + critical_current_gauss[ind_not_zero], + self.dict_Gauss_pt["J_critical"][ind_not_zero], ) + # Evaluate the equivalent electric resistance in Ohm. + self.dict_Gauss_pt["electric_resistance"][ + ind_not_zero] = self.parallel_electric_resistance( + conductor, + [ + "electrical_resistivity_superconductor", + "electrical_resistivity_stabilizer", + ],["sc","stab"], + ind_not_zero, + ) - # Get index of the normal region where all the current is - # carried by the stabilizer, to avoid division by 0 in - # evaluation of superconducting electrical resistivity with the - # power law. - ind_stab_gauss = np.nonzero( - ( - stab_current_gauss / self.dict_Gauss_pt["op_current"][ind_not_zero[ind_normal_gauss]] - > 0.999999 - ) - | (sc_current_gauss < 1.0) - )[0] - - ## CURRENT CARRIED BY THE STABILIZER ## - # Check if np array ind_stab_gauss is not empty. - if ind_stab_gauss.any(): - # Get the index of location of current sharing region, if - # any. - ind_sh_gauss = np.nonzero( - ( - stab_current_gauss - / self.dict_Gauss_pt["op_current"][ind_not_zero[ind_normal_gauss]] - <= 0.999999 - ) - | (sc_current_gauss >= 1.0) - )[0] - - # Check if nparray ind_sh_gauss is not empty. - if ind_sh_gauss.any(): - # ind_sh_gauss is not empty. - # Get final index of the location of the current - # sharing zone, keeping into account that - # sc_current_gauss is already filtered on - # ind_not_zero[ind_normal_gauss]. - ind_shf_gauss = {1:ind_sh_gauss,2:ind_not_zero[ind_normal_gauss[ind_sh_gauss]]} - else: - # ind_sh_gauss is empty. - # Get final index of the location of the current - # sharing zone, keeping into account that - # sc_current_gauss is already filtered on - # ind_not_zero[ind_normal_gauss] and that ind_sh_gauss - # is empty. - ind_shf_gauss = {1:ind_sh_gauss,2:ind_not_zero[ind_normal_gauss]} - - # Evaluate electic resistance in normal region (stabilizer - # only). - self.dict_Gauss_pt["electric_resistance"][ - ind_shf_gauss[2] - ] = self.electric_resistance( - conductor, "electrical_resistivity_stabilizer", "stab", ind_shf_gauss[2] - ) - else: - # Get final index of the location of the current sharing - # zone, keeping into account that sc_current_gauss is - # already filtered on ind_not_zero[ind_normal_gauss]. In - # this case the current is shared between the - # superconductor and the stabilizer, so we need to exploit - # all the index in array ind_normal_gauss for the array - # sc_current_gauss. - ind_shf_gauss = {1:np.nonzero(ind_normal_gauss>=0)[0],2:ind_not_zero[ind_normal_gauss]} - - ## CURRENT SHARED BY THE SUPERCONDUCTOR AND THE STABILIZER ## - if ind_shf_gauss[1].any(): - # Evaluate the electrical resistivity of the superconductor - # according to the power low in Gauss point in Ohm*m. - self.dict_Gauss_pt["electrical_resistivity_superconductor"][ - ind_shf_gauss[2] - ] = self.superconductor_power_law( - sc_current_gauss[ind_shf_gauss[1]], - critical_current_gauss[ind_shf_gauss[2]], - self.dict_Gauss_pt["J_critical"][ind_shf_gauss[2]], - ) - - # Evaluate the equivalent electric resistance in Ohm. - self.dict_Gauss_pt["electric_resistance"][ - ind_shf_gauss[2] - ] = self.parallel_electric_resistance( - conductor, - [ - "electrical_resistivity_superconductor", - "electrical_resistivity_stabilizer", - ],["sc","stab"], - ind_shf_gauss[2], - ) - - # Compute voltage along stabilizer. - v_stab = self.dict_Gauss_pt["electrical_resistivity_stabilizer"][ - ind_shf_gauss[2] - ] * stab_current_gauss[ind_shf_gauss[1]] / self.cross_section["stab"] - # Compute voltage along superconductor - v_sc = self.inputs["E0"] * (sc_current_gauss[ind_shf_gauss[1]] / critical_current_gauss[ind_shf_gauss[2]]) ** self.inputs["nn"] - # Check that the voltage along stabilizer is equal to the - # voltage along superconductor (i.e, check the reliability - # of the current divider). - if all(np.isclose(v_stab,v_sc)) == False: - raise ValueError(f"Voltage difference along superconductor and stabilizer must be the same.") - - return self.dict_Gauss_pt["electric_resistance"] + # Compute voltage along stabilizer. + v_stab = self.dict_Gauss_pt["electrical_resistivity_stabilizer"][ + ind_not_zero + ] * stab_current_gauss / self.cross_section["stab"] + # Compute voltage along superconductor + v_sc = self.inputs["E0"] * (sc_current_gauss / critical_current_gauss[ind_not_zero]) ** self.inputs["nn"] + # Check that the voltage along stabilizer is equal to the + # voltage along superconductor (i.e, check the reliability + # of the current divider). + if all(np.isclose(v_stab,v_sc)) == False: + raise ValueError(f"Voltage difference along superconductor and stabilizer must be the same.") + + return self.dict_Gauss_pt["electric_resistance"] \ No newline at end of file From 20722d604a14fb6d968ace349430c28950455486 Mon Sep 17 00:00:00 2001 From: Mortara Alessandro Date: Mon, 12 Feb 2024 10:57:24 +0100 Subject: [PATCH 06/26] refactor(indices): fun get_electric_resistance Improved the readability of the code and increased its simplicity by removing the multiple indices and distinguishing only two cases: that in which the critical current is =0 and that in which the critical current is other than 0. In the first case, the equivalent resistance will be set equal to the equivalent resistance of the stabilizer part. In the second case, the current divider will be calculated to derive the current on the superconductor, from it with the power law the resistivity of the superconductor will be estimated, and then the equivalent resistance of the stabilizer and superconductor together will be calculated. fun: get_electric_resistance modified: strand_mixed_component.py --- source_code/strand_mixed_component.py | 173 +++++--------------------- 1 file changed, 32 insertions(+), 141 deletions(-) diff --git a/source_code/strand_mixed_component.py b/source_code/strand_mixed_component.py index 76aeb11..350908b 100644 --- a/source_code/strand_mixed_component.py +++ b/source_code/strand_mixed_component.py @@ -797,149 +797,40 @@ def get_electric_resistance(self, conductor: object) -> np.ndarray: # still identify a normal region. if ind_not_zero.any(): - # Get index that correspond to superconducting regime. - ind_sc_gauss = np.nonzero( - self.dict_Gauss_pt["op_current_sc"][ind_not_zero] / critical_current_gauss[ind_not_zero] < 0.95 - )[0] - # Get index that correspond to the normal regime. - ind_normal_gauss = np.nonzero( - self.dict_Gauss_pt["op_current_sc"][ind_not_zero] / critical_current_gauss[ind_not_zero] >= 0.95 - )[0] - - ## SUPERCONDUCTING REGIME ## - - # Check if np array ind_sc_gauss is not empty. - if ind_sc_gauss.any(): - # Current in superconducting regime is the carried by - # superconducting material only. - self.dict_Gauss_pt["op_current"][ind_not_zero[ind_sc_gauss]] = self.dict_Gauss_pt[ - "op_current_sc" - ][ind_not_zero[ind_sc_gauss]] - - # Compute superconducting electrical resistivity only in index - # for which the superconducting regime is guaranteed, using the - # power low. - self.dict_Gauss_pt["electrical_resistivity_superconductor"][ind_not_zero[ind_sc_gauss]] = self.superconductor_power_law( - self.dict_Gauss_pt["op_current"][ind_not_zero[ind_sc_gauss]], - critical_current_gauss[ind_not_zero[ind_sc_gauss]], - self.dict_Gauss_pt["J_critical"][ind_not_zero[ind_sc_gauss]] + sc_current_gauss, stab_current_gauss = self.solve_current_divider( + self.dict_Gauss_pt["electrical_resistivity_stabilizer"][ind_not_zero], + critical_current_gauss[ind_not_zero], + self.dict_Gauss_pt["op_current"][ind_not_zero] ) - - # Evaluate electic resistance in superconducting region - # (superconductor only). - self.dict_Gauss_pt["electric_resistance"][ind_not_zero[ - ind_sc_gauss - ]] = self.electric_resistance( - conductor, "electrical_resistivity_superconductor", "sc", ind_not_zero[ind_sc_gauss] - ) - - ## NORMAL REGIME ## - - # Check if np array ind_normal_gauss is not empty. - if ind_normal_gauss.any(): - - # Evaluate how the current is distributed solving the current - # divider problem in Gauss point. - sc_current_gauss, stab_current_gauss = self.solve_current_divider( - self.dict_Gauss_pt["electrical_resistivity_stabilizer"][ind_not_zero[ind_normal_gauss]], - critical_current_gauss[ind_not_zero[ind_normal_gauss]], - self.dict_Gauss_pt["op_current"][ind_not_zero[ind_normal_gauss]] + + self.dict_Gauss_pt["electrical_resistivity_superconductor"][ + ind_not_zero + ] = self.superconductor_power_law( + sc_current_gauss, + critical_current_gauss[ind_not_zero], + self.dict_Gauss_pt["J_critical"][ind_not_zero], ) + # Evaluate the equivalent electric resistance in Ohm. + self.dict_Gauss_pt["electric_resistance"][ + ind_not_zero] = self.parallel_electric_resistance( + conductor, + [ + "electrical_resistivity_superconductor", + "electrical_resistivity_stabilizer", + ],["sc","stab"], + ind_not_zero, + ) - # Get index of the normal region where all the current is - # carried by the stabilizer, to avoid division by 0 in - # evaluation of superconducting electrical resistivity with the - # power law. - ind_stab_gauss = np.nonzero( - ( - stab_current_gauss / self.dict_Gauss_pt["op_current"][ind_not_zero[ind_normal_gauss]] - > 0.999999 - ) - | (sc_current_gauss < 1.0) - )[0] - - ## CURRENT CARRIED BY THE STABILIZER ## - # Check if np array ind_stab_gauss is not empty. - if ind_stab_gauss.any(): - # Get the index of location of current sharing region, if - # any. - ind_sh_gauss = np.nonzero( - ( - stab_current_gauss - / self.dict_Gauss_pt["op_current"][ind_not_zero[ind_normal_gauss]] - <= 0.999999 - ) - | (sc_current_gauss >= 1.0) - )[0] - - # Check if nparray ind_sh_gauss is not empty. - if ind_sh_gauss.any(): - # ind_sh_gauss is not empty. - # Get final index of the location of the current - # sharing zone, keeping into account that - # sc_current_gauss is already filtered on - # ind_not_zero[ind_normal_gauss]. - ind_shf_gauss = {1:ind_sh_gauss,2:ind_not_zero[ind_normal_gauss[ind_sh_gauss]]} - else: - # ind_sh_gauss is empty. - # Get final index of the location of the current - # sharing zone, keeping into account that - # sc_current_gauss is already filtered on - # ind_not_zero[ind_normal_gauss] and that ind_sh_gauss - # is empty. - ind_shf_gauss = {1:ind_sh_gauss,2:ind_not_zero[ind_normal_gauss]} - - # Evaluate electic resistance in normal region (stabilizer - # only). - self.dict_Gauss_pt["electric_resistance"][ - ind_shf_gauss[2] - ] = self.electric_resistance( - conductor, "electrical_resistivity_stabilizer", "stab", ind_shf_gauss[2] - ) - else: - # Get final index of the location of the current sharing - # zone, keeping into account that sc_current_gauss is - # already filtered on ind_not_zero[ind_normal_gauss]. In - # this case the current is shared between the - # superconductor and the stabilizer, so we need to exploit - # all the index in array ind_normal_gauss for the array - # sc_current_gauss. - ind_shf_gauss = {1:np.nonzero(ind_normal_gauss>=0)[0],2:ind_not_zero[ind_normal_gauss]} - - ## CURRENT SHARED BY THE SUPERCONDUCTOR AND THE STABILIZER ## - if ind_shf_gauss[1].any(): - # Evaluate the electrical resistivity of the superconductor - # according to the power low in Gauss point in Ohm*m. - self.dict_Gauss_pt["electrical_resistivity_superconductor"][ - ind_shf_gauss[2] - ] = self.superconductor_power_law( - sc_current_gauss[ind_shf_gauss[1]], - critical_current_gauss[ind_shf_gauss[2]], - self.dict_Gauss_pt["J_critical"][ind_shf_gauss[2]], - ) - - # Evaluate the equivalent electric resistance in Ohm. - self.dict_Gauss_pt["electric_resistance"][ - ind_shf_gauss[2] - ] = self.parallel_electric_resistance( - conductor, - [ - "electrical_resistivity_superconductor", - "electrical_resistivity_stabilizer", - ],["sc","stab"], - ind_shf_gauss[2], - ) - - # Compute voltage along stabilizer. - v_stab = self.dict_Gauss_pt["electrical_resistivity_stabilizer"][ - ind_shf_gauss[2] - ] * stab_current_gauss[ind_shf_gauss[1]] / self.cross_section["stab"] - # Compute voltage along superconductor - v_sc = self.inputs["E0"] * (sc_current_gauss[ind_shf_gauss[1]] / critical_current_gauss[ind_shf_gauss[2]]) ** self.inputs["nn"] - # Check that the voltage along stabilizer is equal to the - # voltage along superconductor (i.e, check the reliability - # of the current divider). - if all(np.isclose(v_stab,v_sc)) == False: - raise ValueError(f"Voltage difference along superconductor and stabilizer must be the same.") + # Compute voltage along stabilizer. + v_stab = self.dict_Gauss_pt["electrical_resistivity_stabilizer"][ + ind_not_zero + ] * stab_current_gauss / self.cross_section["stab"] + # Compute voltage along superconductor + v_sc = self.inputs["E0"] * (sc_current_gauss / critical_current_gauss[ind_not_zero]) ** self.inputs["nn"] + # Check that the voltage along stabilizer is equal to the + # voltage along superconductor (i.e, check the reliability + # of the current divider). + if all(np.isclose(v_stab,v_sc)) == False: + raise ValueError(f"Voltage difference along superconductor and stabilizer must be the same.") return self.dict_Gauss_pt["electric_resistance"] \ No newline at end of file From fbad580ad0b8aa7823e5f1ab40f0b7563518f78c Mon Sep 17 00:00:00 2001 From: Mortara Alessandro Date: Mon, 12 Feb 2024 19:21:44 +0100 Subject: [PATCH 07/26] chore: restored the inductance methods The inductance methods were commented and a constant value was imposed, to confront the results with external references. All the calculation of both the self and the mutual inductance are now evaluated by the code. modified: conductor.py --- source_code/conductor.py | 229 +++++++++++++++++++-------------------- 1 file changed, 113 insertions(+), 116 deletions(-) diff --git a/source_code/conductor.py b/source_code/conductor.py index 0b1e8eb..a356890 100644 --- a/source_code/conductor.py +++ b/source_code/conductor.py @@ -3776,92 +3776,91 @@ def __mutual_inductance( """ jj = np.r_[ii + 1 : self.total_elements_current_carriers] - # ll = lmod[jj] - # mm = lmod[ii] - # len_jj = self.total_elements_current_carriers - (ii + 1) - # rr = dict( - # end_end=np.zeros(len_jj), - # end_start=np.zeros(len_jj), - # start_start=np.zeros(len_jj), - # start_end=np.zeros(len_jj), - # ) + ll = lmod[jj] + mm = lmod[ii] + len_jj = self.total_elements_current_carriers - (ii + 1) + rr = dict( + end_end=np.zeros(len_jj), + end_start=np.zeros(len_jj), + start_start=np.zeros(len_jj), + start_end=np.zeros(len_jj), + ) - # for key in rr.keys(): - # rr[key] = self.__vertex_to_vertex_distance(key, ii, jj) - # # End for key + for key in rr.keys(): + rr[key] = self.__vertex_to_vertex_distance(key, ii, jj) + # End for key - # # Additional parameters - # alpha2 = ( - # rr["start_end"] ** 2 - # - rr["start_start"] ** 2 - # + rr["end_start"] ** 2 - # - rr["end_end"] ** 2 - # ) + # Additional parameters + alpha2 = ( + rr["start_end"] ** 2 + - rr["start_start"] ** 2 + + rr["end_start"] ** 2 + - rr["end_end"] ** 2 + ) - # cos_eps = np.minimum(np.maximum(alpha2 / (2 * ll * mm), -1.0), 1.0) - # sin_eps = np.sin(np.arccos(cos_eps)) - - # dd = 4 * ll ** 2 * mm ** 2 - alpha2 ** 2 - # mu = ( - # ll - # * ( - # 2 * mm ** 2 * (rr["end_start"] ** 2 - rr["start_start"] ** 2 - ll ** 2) - # + alpha2 * (rr["start_end"] ** 2 - rr["start_start"] ** 2 - mm ** 2) - # ) - # / dd - # ) - # nu = ( - # mm - # * ( - # 2 * ll ** 2 * (rr["start_end"] ** 2 - rr["start_start"] ** 2 - mm ** 2) - # + alpha2 * (rr["end_start"] ** 2 - rr["start_start"] ** 2 - ll ** 2) - # ) - # / dd - # ) - # d2 = rr["start_start"] ** 2 - mu ** 2 - nu ** 2 + 2 * mu * nu * cos_eps - - # # avoid rounding for segments in a plane - # d2[d2 < abstol ** 2] = 0 - # d0 = np.sqrt(d2) - - # # solid angles - # omega = ( - # np.arctan( - # (d2 * cos_eps + (mu + ll) * (nu + mm) * sin_eps ** 2) - # / (d0 * rr["end_end"] * sin_eps) - # ) - # - np.arctan( - # (d2 * cos_eps + (mu + ll) * nu * sin_eps ** 2) - # / (d0 * rr["end_start"] * sin_eps) - # ) - # + np.arctan( - # (d2 * cos_eps + mu * nu * sin_eps ** 2) - # / (d0 * rr["start_start"] * sin_eps) - # ) - # - np.arctan( - # (d2 * cos_eps + mu * (nu + mm) * sin_eps ** 2) - # / (d0 * rr["start_end"] * sin_eps) - # ) - # ) - # omega[d0 == 0.0] = 0.0 + cos_eps = np.minimum(np.maximum(alpha2 / (2 * ll * mm), -1.0), 1.0) + sin_eps = np.sin(np.arccos(cos_eps)) - # # contribution - # pp = np.zeros((len_jj, 5), dtype=float) - # pp[:, 0] = (ll + mu) * np.arctanh(mm / (rr["end_end"] + rr["end_start"])) - # pp[:, 1] = -nu * np.arctanh(ll / (rr["end_start"] + rr["start_start"])) - # pp[:, 2] = (mm + nu) * np.arctanh(ll / (rr["end_end"] + rr["start_end"])) - # pp[:, 3] = -mu * np.arctanh(mm / (rr["start_start"] + rr["start_end"])) - # pp[:, 4] = d0 * omega / sin_eps + dd = 4 * ll ** 2 * mm ** 2 - alpha2 ** 2 + mu = ( + ll + * ( + 2 * mm ** 2 * (rr["end_start"] ** 2 - rr["start_start"] ** 2 - ll ** 2) + + alpha2 * (rr["start_end"] ** 2 - rr["start_start"] ** 2 - mm ** 2) + ) + / dd + ) + nu = ( + mm + * ( + 2 * ll ** 2 * (rr["start_end"] ** 2 - rr["start_start"] ** 2 - mm ** 2) + + alpha2 * (rr["end_start"] ** 2 - rr["start_start"] ** 2 - ll ** 2) + ) + / dd + ) + d2 = rr["start_start"] ** 2 - mu ** 2 - nu ** 2 + 2 * mu * nu * cos_eps - # # filter odd cases (e.g. consecutive segments) - # pp[np.isnan(pp)] = 0.0 - # pp[np.isinf(pp)] = 0.0 + # avoid rounding for segments in a plane + d2[d2 < abstol ** 2] = 0 + d0 = np.sqrt(d2) + + # solid angles + omega = ( + np.arctan( + (d2 * cos_eps + (mu + ll) * (nu + mm) * sin_eps ** 2) + / (d0 * rr["end_end"] * sin_eps) + ) + - np.arctan( + (d2 * cos_eps + (mu + ll) * nu * sin_eps ** 2) + / (d0 * rr["end_start"] * sin_eps) + ) + + np.arctan( + (d2 * cos_eps + mu * nu * sin_eps ** 2) + / (d0 * rr["start_start"] * sin_eps) + ) + - np.arctan( + (d2 * cos_eps + mu * (nu + mm) * sin_eps ** 2) + / (d0 * rr["start_end"] * sin_eps) + ) + ) + omega[d0 == 0.0] = 0.0 + + # contribution + pp = np.zeros((len_jj, 5), dtype=float) + pp[:, 0] = (ll + mu) * np.arctanh(mm / (rr["end_end"] + rr["end_start"])) + pp[:, 1] = -nu * np.arctanh(ll / (rr["end_start"] + rr["start_start"])) + pp[:, 2] = (mm + nu) * np.arctanh(ll / (rr["end_end"] + rr["start_end"])) + pp[:, 3] = -mu * np.arctanh(mm / (rr["start_start"] + rr["start_end"])) + pp[:, 4] = d0 * omega / sin_eps + + # filter odd cases (e.g. consecutive segments) + pp[np.isnan(pp)] = 0.0 + pp[np.isinf(pp)] = 0.0 # Mutual inductances matrix[ii, jj] = ( - 5.0e-8 # Mutual inductance imposed by Zappatore input file - # 2 * cos_eps * (pp[:, 0] + pp[:, 1] + pp[:, 2] + pp[:, 3]) - # - cos_eps * pp[:, 4] + 2 * cos_eps * (pp[:, 0] + pp[:, 1] + pp[:, 2] + pp[:, 3]) + - cos_eps * pp[:, 4] ) return matrix @@ -3921,24 +3920,23 @@ def __self_inductance_mode1(self, lmod: np.ndarray) -> np.ndarray: for ii, obj in enumerate(self.inventory["StrandComponent"].collection): self_inductance[ii :: self.inventory["StrandComponent"].number] = ( - 1.0e-7 - # 2 - # * lmod[ii :: self.inventory["StrandComponent"].number] - # * ( - # np.arcsinh( - # lmod[ii :: self.inventory["StrandComponent"].number] - # / obj.radius - # ) - # - np.sqrt( - # 1.0 - # + ( - # obj.radius - # / lmod[ii :: self.inventory["StrandComponent"].number] - # ) - # ** 2 - # ) - # + obj.radius / lmod[ii :: self.inventory["StrandComponent"].number] - # ) + 2 + * lmod[ii :: self.inventory["StrandComponent"].number] + * ( + np.arcsinh( + lmod[ii :: self.inventory["StrandComponent"].number] + / obj.radius + ) + - np.sqrt( + 1.0 + + ( + obj.radius + / lmod[ii :: self.inventory["StrandComponent"].number] + ) + ** 2 + ) + + obj.radius / lmod[ii :: self.inventory["StrandComponent"].number] + ) ) return self_inductance @@ -3955,26 +3953,25 @@ def __self_inductance_mode2(self, lmod: np.ndarray) -> np.ndarray: for ii, obj in enumerate(self.inventory["StrandComponent"].collection): self_inductance[ii :: self.inventory["StrandComponent"].number] = ( - 1.0e-7 -# 2 * ( -# lmod[ii :: self.inventory["StrandComponent"].number] -# * np.log( -# ( -# lmod[ii :: self.inventory["StrandComponent"].number] -# + np.sqrt( -# lmod[ii :: self.inventory["StrandComponent"].number] ** 2 -# + obj.radius ** 2 -# ) -# ) -# / obj.radius -# ) -# - np.sqrt( -# lmod[ii :: self.inventory["StrandComponent"].number] ** 2 -# + obj.radius ** 2 -# ) -# + lmod[ii :: self.inventory["StrandComponent"].number] / 4 -# + obj.radius - ) + 2 * ( + lmod[ii :: self.inventory["StrandComponent"].number] + * np.log( + ( + lmod[ii :: self.inventory["StrandComponent"].number] + + np.sqrt( + lmod[ii :: self.inventory["StrandComponent"].number] ** 2 + + obj.radius ** 2 + ) + ) + / obj.radius + ) + - np.sqrt( + lmod[ii :: self.inventory["StrandComponent"].number] ** 2 + + obj.radius ** 2 + ) + + lmod[ii :: self.inventory["StrandComponent"].number] / 4 + + obj.radius + )) return self_inductance From 51cfa7e1d4baa0b247ebb990c03ff703fe61f66d Mon Sep 17 00:00:00 2001 From: Mortara Alessandro Date: Mon, 12 Feb 2024 19:32:10 +0100 Subject: [PATCH 08/26] chore: new flags for inductances evaluation Two new flags were added: * SELF_INDUCTANCE_MODE_0 = 0 is used to assign a constant value for the self inductance readed from input file conductor_definition.xlsx * CONSTANT_INDUCTANCE = 0 is used to assign a constant value for the mutual inductance readed from input file conductor_definition.xlsx The flags related to mutual inductance were redefined as: * ANALYTICAL_INDUCTANCE = 1 Flag to evaluate inductance analytically * APPROXIMATE_INDUCTANCE = 2 Flag to evaluate inductance using an approximation modified: conductor_flags.py --- source_code/conductor_flags.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/source_code/conductor_flags.py b/source_code/conductor_flags.py index 7c61647..4d87e10 100644 --- a/source_code/conductor_flags.py +++ b/source_code/conductor_flags.py @@ -16,15 +16,19 @@ # Electric conductance is not defined per unit length ELECTRIC_CONDUCTANCE_NOT_UNIT_LENGTH = 2 +# Self inductance readed from input file conductor_definition.xlsx +SELF_INDUCTANCE_MODE_0 = 0 # Analytical self inductance is evaluated according to mode 1. SELF_INDUCTANCE_MODE_1 = 1 # Analytical self inductance is evaluated according to mode 2. SELF_INDUCTANCE_MODE_2 = 2 +# Constant mutual inductance readed from input file conductor_definition.xlsx +CONSTANT_INDUCTANCE = 0 # Flag to evaluate inductance analytically -ANALYTICAL_INDUCTANCE = 0 +ANALYTICAL_INDUCTANCE = 1 # Flag to evaluate inductance using an approximation -APPROXIMATE_INDUCTANCE = 1 +APPROXIMATE_INDUCTANCE = 2 # Flag to solve the electric problem in steady state conditions. STATIC_ELECTRIC_SOLVER = 0 From c01a0d1187fc318c9a344c5643ff2c05573fda6c Mon Sep 17 00:00:00 2001 From: Mortara Alessandro Date: Tue, 13 Feb 2024 10:49:59 +0100 Subject: [PATCH 09/26] Revert "chore: restored the inductance methods" This reverts commit fbad580ad0b8aa7823e5f1ab40f0b7563518f78c. --- source_code/conductor.py | 229 ++++++++++++++++++++------------------- 1 file changed, 116 insertions(+), 113 deletions(-) diff --git a/source_code/conductor.py b/source_code/conductor.py index a356890..0b1e8eb 100644 --- a/source_code/conductor.py +++ b/source_code/conductor.py @@ -3776,91 +3776,92 @@ def __mutual_inductance( """ jj = np.r_[ii + 1 : self.total_elements_current_carriers] - ll = lmod[jj] - mm = lmod[ii] - len_jj = self.total_elements_current_carriers - (ii + 1) - rr = dict( - end_end=np.zeros(len_jj), - end_start=np.zeros(len_jj), - start_start=np.zeros(len_jj), - start_end=np.zeros(len_jj), - ) - - for key in rr.keys(): - rr[key] = self.__vertex_to_vertex_distance(key, ii, jj) - # End for key - - # Additional parameters - alpha2 = ( - rr["start_end"] ** 2 - - rr["start_start"] ** 2 - + rr["end_start"] ** 2 - - rr["end_end"] ** 2 - ) + # ll = lmod[jj] + # mm = lmod[ii] + # len_jj = self.total_elements_current_carriers - (ii + 1) + # rr = dict( + # end_end=np.zeros(len_jj), + # end_start=np.zeros(len_jj), + # start_start=np.zeros(len_jj), + # start_end=np.zeros(len_jj), + # ) - cos_eps = np.minimum(np.maximum(alpha2 / (2 * ll * mm), -1.0), 1.0) - sin_eps = np.sin(np.arccos(cos_eps)) + # for key in rr.keys(): + # rr[key] = self.__vertex_to_vertex_distance(key, ii, jj) + # # End for key - dd = 4 * ll ** 2 * mm ** 2 - alpha2 ** 2 - mu = ( - ll - * ( - 2 * mm ** 2 * (rr["end_start"] ** 2 - rr["start_start"] ** 2 - ll ** 2) - + alpha2 * (rr["start_end"] ** 2 - rr["start_start"] ** 2 - mm ** 2) - ) - / dd - ) - nu = ( - mm - * ( - 2 * ll ** 2 * (rr["start_end"] ** 2 - rr["start_start"] ** 2 - mm ** 2) - + alpha2 * (rr["end_start"] ** 2 - rr["start_start"] ** 2 - ll ** 2) - ) - / dd - ) - d2 = rr["start_start"] ** 2 - mu ** 2 - nu ** 2 + 2 * mu * nu * cos_eps - - # avoid rounding for segments in a plane - d2[d2 < abstol ** 2] = 0 - d0 = np.sqrt(d2) + # # Additional parameters + # alpha2 = ( + # rr["start_end"] ** 2 + # - rr["start_start"] ** 2 + # + rr["end_start"] ** 2 + # - rr["end_end"] ** 2 + # ) - # solid angles - omega = ( - np.arctan( - (d2 * cos_eps + (mu + ll) * (nu + mm) * sin_eps ** 2) - / (d0 * rr["end_end"] * sin_eps) - ) - - np.arctan( - (d2 * cos_eps + (mu + ll) * nu * sin_eps ** 2) - / (d0 * rr["end_start"] * sin_eps) - ) - + np.arctan( - (d2 * cos_eps + mu * nu * sin_eps ** 2) - / (d0 * rr["start_start"] * sin_eps) - ) - - np.arctan( - (d2 * cos_eps + mu * (nu + mm) * sin_eps ** 2) - / (d0 * rr["start_end"] * sin_eps) - ) - ) - omega[d0 == 0.0] = 0.0 + # cos_eps = np.minimum(np.maximum(alpha2 / (2 * ll * mm), -1.0), 1.0) + # sin_eps = np.sin(np.arccos(cos_eps)) + + # dd = 4 * ll ** 2 * mm ** 2 - alpha2 ** 2 + # mu = ( + # ll + # * ( + # 2 * mm ** 2 * (rr["end_start"] ** 2 - rr["start_start"] ** 2 - ll ** 2) + # + alpha2 * (rr["start_end"] ** 2 - rr["start_start"] ** 2 - mm ** 2) + # ) + # / dd + # ) + # nu = ( + # mm + # * ( + # 2 * ll ** 2 * (rr["start_end"] ** 2 - rr["start_start"] ** 2 - mm ** 2) + # + alpha2 * (rr["end_start"] ** 2 - rr["start_start"] ** 2 - ll ** 2) + # ) + # / dd + # ) + # d2 = rr["start_start"] ** 2 - mu ** 2 - nu ** 2 + 2 * mu * nu * cos_eps + + # # avoid rounding for segments in a plane + # d2[d2 < abstol ** 2] = 0 + # d0 = np.sqrt(d2) + + # # solid angles + # omega = ( + # np.arctan( + # (d2 * cos_eps + (mu + ll) * (nu + mm) * sin_eps ** 2) + # / (d0 * rr["end_end"] * sin_eps) + # ) + # - np.arctan( + # (d2 * cos_eps + (mu + ll) * nu * sin_eps ** 2) + # / (d0 * rr["end_start"] * sin_eps) + # ) + # + np.arctan( + # (d2 * cos_eps + mu * nu * sin_eps ** 2) + # / (d0 * rr["start_start"] * sin_eps) + # ) + # - np.arctan( + # (d2 * cos_eps + mu * (nu + mm) * sin_eps ** 2) + # / (d0 * rr["start_end"] * sin_eps) + # ) + # ) + # omega[d0 == 0.0] = 0.0 - # contribution - pp = np.zeros((len_jj, 5), dtype=float) - pp[:, 0] = (ll + mu) * np.arctanh(mm / (rr["end_end"] + rr["end_start"])) - pp[:, 1] = -nu * np.arctanh(ll / (rr["end_start"] + rr["start_start"])) - pp[:, 2] = (mm + nu) * np.arctanh(ll / (rr["end_end"] + rr["start_end"])) - pp[:, 3] = -mu * np.arctanh(mm / (rr["start_start"] + rr["start_end"])) - pp[:, 4] = d0 * omega / sin_eps + # # contribution + # pp = np.zeros((len_jj, 5), dtype=float) + # pp[:, 0] = (ll + mu) * np.arctanh(mm / (rr["end_end"] + rr["end_start"])) + # pp[:, 1] = -nu * np.arctanh(ll / (rr["end_start"] + rr["start_start"])) + # pp[:, 2] = (mm + nu) * np.arctanh(ll / (rr["end_end"] + rr["start_end"])) + # pp[:, 3] = -mu * np.arctanh(mm / (rr["start_start"] + rr["start_end"])) + # pp[:, 4] = d0 * omega / sin_eps - # filter odd cases (e.g. consecutive segments) - pp[np.isnan(pp)] = 0.0 - pp[np.isinf(pp)] = 0.0 + # # filter odd cases (e.g. consecutive segments) + # pp[np.isnan(pp)] = 0.0 + # pp[np.isinf(pp)] = 0.0 # Mutual inductances matrix[ii, jj] = ( - 2 * cos_eps * (pp[:, 0] + pp[:, 1] + pp[:, 2] + pp[:, 3]) - - cos_eps * pp[:, 4] + 5.0e-8 # Mutual inductance imposed by Zappatore input file + # 2 * cos_eps * (pp[:, 0] + pp[:, 1] + pp[:, 2] + pp[:, 3]) + # - cos_eps * pp[:, 4] ) return matrix @@ -3920,23 +3921,24 @@ def __self_inductance_mode1(self, lmod: np.ndarray) -> np.ndarray: for ii, obj in enumerate(self.inventory["StrandComponent"].collection): self_inductance[ii :: self.inventory["StrandComponent"].number] = ( - 2 - * lmod[ii :: self.inventory["StrandComponent"].number] - * ( - np.arcsinh( - lmod[ii :: self.inventory["StrandComponent"].number] - / obj.radius - ) - - np.sqrt( - 1.0 - + ( - obj.radius - / lmod[ii :: self.inventory["StrandComponent"].number] - ) - ** 2 - ) - + obj.radius / lmod[ii :: self.inventory["StrandComponent"].number] - ) + 1.0e-7 + # 2 + # * lmod[ii :: self.inventory["StrandComponent"].number] + # * ( + # np.arcsinh( + # lmod[ii :: self.inventory["StrandComponent"].number] + # / obj.radius + # ) + # - np.sqrt( + # 1.0 + # + ( + # obj.radius + # / lmod[ii :: self.inventory["StrandComponent"].number] + # ) + # ** 2 + # ) + # + obj.radius / lmod[ii :: self.inventory["StrandComponent"].number] + # ) ) return self_inductance @@ -3953,25 +3955,26 @@ def __self_inductance_mode2(self, lmod: np.ndarray) -> np.ndarray: for ii, obj in enumerate(self.inventory["StrandComponent"].collection): self_inductance[ii :: self.inventory["StrandComponent"].number] = ( - 2 * ( - lmod[ii :: self.inventory["StrandComponent"].number] - * np.log( - ( - lmod[ii :: self.inventory["StrandComponent"].number] - + np.sqrt( - lmod[ii :: self.inventory["StrandComponent"].number] ** 2 - + obj.radius ** 2 - ) - ) - / obj.radius - ) - - np.sqrt( - lmod[ii :: self.inventory["StrandComponent"].number] ** 2 - + obj.radius ** 2 - ) - + lmod[ii :: self.inventory["StrandComponent"].number] / 4 - + obj.radius - )) + 1.0e-7 +# 2 * ( +# lmod[ii :: self.inventory["StrandComponent"].number] +# * np.log( +# ( +# lmod[ii :: self.inventory["StrandComponent"].number] +# + np.sqrt( +# lmod[ii :: self.inventory["StrandComponent"].number] ** 2 +# + obj.radius ** 2 +# ) +# ) +# / obj.radius +# ) +# - np.sqrt( +# lmod[ii :: self.inventory["StrandComponent"].number] ** 2 +# + obj.radius ** 2 +# ) +# + lmod[ii :: self.inventory["StrandComponent"].number] / 4 +# + obj.radius + ) return self_inductance From b1aad2fb7e372d35eebc0b5d50b9766c0fd609e2 Mon Sep 17 00:00:00 2001 From: Mortara Alessandro Date: Tue, 13 Feb 2024 12:57:33 +0100 Subject: [PATCH 10/26] feat!: new methods in class Conductor News privates methods to deal with the inductance evaluation. *__constant_inductance: private method that assigns a constant value to the mutual inductance as defined by the user in variable MUTUAL_INDUCTANCE, in sheet CONDUCTOR_operation of the input file conductor_definition.xlsx. *__constant_self_inductance_evaluation: private method that assigns a constant value to the self inductance that is defined by the user in variable SELF_INDUCTANCE, in sheet CONDICTOR_operation in the input file conductor_definition.xlsx Class: Conductor modified: conductor.py BREAKING CHANGES These methods are called respectively with the new values of flags INDUCTANCE_MODE, SELF_INDUCTANCE_MODE, as follows: INDUCTANCE_MODE = 0 -> __constant_inductance SELF_INDUCTANCE_MODE = 0 -> __constant_self_inductance_evaluation --- source_code/conductor.py | 71 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) diff --git a/source_code/conductor.py b/source_code/conductor.py index 0b1e8eb..3fee5fc 100644 --- a/source_code/conductor.py +++ b/source_code/conductor.py @@ -3688,6 +3688,62 @@ def build_electric_known_term_vector(self): self.total_elements_current_carriers : ] = self.dict_node_pt["op_current"] + # CONSTANT INDUCTANCE + def __constant_inductance(self, mode: int): + """Private method that assigns a constant value to the mutual inductance as defined by the user in sheet CONDUCTOR_operation of the input file conductor_definition.xlsx + + Args: + mode (int): flag to select the equation for the analytical evaluation of self inductance. 0:constan value from sheet CONDUCTOR_operation of the input file conductor_definition.xlsx; 1: from method __self_inductance_mode1; 2: from method __self_inductance_mode2. + """ + + lmod = ( + ( + ( + self.nodal_coordinates.iloc[ + self.connectivity_matrix.loc[ + "StrandComponent", + "end", + ], + :, + ] + - self.nodal_coordinates.iloc[ + self.connectivity_matrix.loc[ + "StrandComponent", + "start", + ], + :, + ] + ) + ** 2 + ) + .sum(axis=1) + .apply(np.sqrt) + ) + mutual_inductance = self.operations["MUTUAL_INDUCTANCE"] * np.ones(self.inductance_matrix.shape) + + # The principal diagonal is set to 0 + for ii in range(mutual_inductance.shape[0]): + mutual_inductance[ii, ii] = 0 + + self_inductance_switch = { + SELF_INDUCTANCE_MODE_0: self.__constant_self_inductance_evaluation, + SELF_INDUCTANCE_MODE_1: self.__self_inductance_mode1, + SELF_INDUCTANCE_MODE_2: self.__self_inductance_mode2, + } + self_inductance = self_inductance_switch[mode](lmod) + + # Evaluate internal inductance + internal_inductance = lmod.to_numpy() / 2.0 + + self.inductance_matrix = ( + constants.mu_0 + / (4.0 * constants.pi) + * ( + np.diag(self_inductance + internal_inductance) + + mutual_inductance + + mutual_inductance.T + ) + ) # START: INDUCTANCE ANALYTICAL EVALUATION def __inductance_analytical_calculation(self, mode: int = 2): @@ -3908,6 +3964,21 @@ def __vertex_to_vertex_distance( .apply(np.sqrt) ) + # CONSTANT SELF INDUCTANCE + def __constant_self_inductance_evaluation(self, lmod: np.array) -> np.ndarray: + """Private method that assigns a constant value to the self inductance that is defined by the user in sheet CONDICTOR_operation in the input file conductor_definition.xlsx + + Args: + lmod (np.ndarray): array with the distance between strand component nodal nodes. + + Returns: + np.ndarray: self inductances. + """ + + self_inductance = np.ones(lmod.shape) * self.operations["SELF_INDUCTANCE"] + + return self_inductance + def __self_inductance_mode1(self, lmod: np.ndarray) -> np.ndarray: """Private method that analytically evaluates self inductances according to mode 1. From 740092a48e8557a6d671ee241437bad5984a9588 Mon Sep 17 00:00:00 2001 From: Mortara Alessandro Date: Tue, 13 Feb 2024 14:33:51 +0100 Subject: [PATCH 11/26] chore: update methods for inductance evaluation Modified functions: *__inductance_analytical_calculation *added check on the new mode (value = 0) *added keyword SELF_INDUCTANCE_MODE_0=0 in self_inductance_switch __inductance_approximate_calculation *added check on the new mode (value = 0) *added keyword SELF_INDUCTANCE_MODE_0=0 in self_inductance_switch __build_electric_mass_matrix *added check on the new mode (value = 0) *added keyword CONSTANT_INDUCTANCE=0 in inductance_switch NEED REFACTORING!!! Class: Conductor modified: conductor.py --- source_code/conductor.py | 33 ++++++++++++++++++++++++--------- 1 file changed, 24 insertions(+), 9 deletions(-) diff --git a/source_code/conductor.py b/source_code/conductor.py index 3fee5fc..440ae31 100644 --- a/source_code/conductor.py +++ b/source_code/conductor.py @@ -14,6 +14,7 @@ # import classes from component_collection import ComponentCollection from conductor_flags import ( + CONSTANT_INDUCTANCE, ANALYTICAL_INDUCTANCE, APPROXIMATE_INDUCTANCE, ELECTRIC_CONDUCTANCE_UNIT_LENGTH, @@ -22,6 +23,7 @@ IOP_CONSTANT, IOP_FROM_EXT_FUNCTION, IOP_FROM_FILE, + SELF_INDUCTANCE_MODE_0, SELF_INDUCTANCE_MODE_1, SELF_INDUCTANCE_MODE_2, STATIC_ELECTRIC_SOLVER, @@ -3753,9 +3755,9 @@ def __inductance_analytical_calculation(self, mode: int = 2): mode (int,optional): flag to select the equation for the analytical evaluation of self inductance. 1: mode 1; 2: mode 2. Defaults to 2. """ - if mode != 1 and mode != 2: + if mode!= SELF_INDUCTANCE_MODE_0 and mode != SELF_INDUCTANCE_MODE_1 and mode != SELF_INDUCTANCE_MODE_2: raise ValueError( - f"{self.identifier}\nArgument 'mode' must be equal to {SELF_INDUCTANCE_MODE_1 = } or to {SELF_INDUCTANCE_MODE_2 = }. Current value {mode = } is not allowed. Please check sheet {self.workbook_sheet_name[2]} in file {self.workbook_name}.\n" + f"{self.identifier}\nArgument 'mode' must be equal to {SELF_INDUCTANCE_MODE_0 = } or to {SELF_INDUCTANCE_MODE_1 = } or to {SELF_INDUCTANCE_MODE_2 = }. Current value {mode = } is not allowed. Please check sheet {self.workbook_sheet_name[2]} in file {self.workbook_name}.\n" ) ABSTOL = 1e-6 lmod = ( @@ -3799,6 +3801,7 @@ def __inductance_analytical_calculation(self, mode: int = 2): # Switch to evalutae self inductance. self_inductance_switch = { + SELF_INDUCTANCE_MODE_0: self.__constant_self_inductance_evaluation, SELF_INDUCTANCE_MODE_1: self.__self_inductance_mode1, SELF_INDUCTANCE_MODE_2: self.__self_inductance_mode2, } @@ -4053,8 +4056,13 @@ def __self_inductance_mode2(self, lmod: np.ndarray) -> np.ndarray: # START: INDUCTANCE APPROXIMATE EVALUATION - def __inductance_approximate_calculation(self): + def __inductance_approximate_calculation(self, mode : int = 2): """Private method that approximate the inductance of the system. For an analytical evaluation of the inductance use private method __inductance_analytical_calculation.""" + + if mode!= SELF_INDUCTANCE_MODE_0 and mode != SELF_INDUCTANCE_MODE_1 and mode != SELF_INDUCTANCE_MODE_2: + raise ValueError( + f"{self.identifier}\nArgument 'mode' must be equal to {SELF_INDUCTANCE_MODE_0 = } or to {SELF_INDUCTANCE_MODE_1 = } or to {SELF_INDUCTANCE_MODE_2 = }. Current value {mode = } is not allowed. Please check sheet {self.workbook_sheet_name[2]} in file {self.workbook_name}.\n" + ) ll = ( self.nodal_coordinates.iloc[ @@ -4088,8 +4096,13 @@ def __inductance_approximate_calculation(self): ll, mutual_inductance ) - # Evaluate self inductance - self_inductance = self.__self_inductance_approximate(lmod) + self_inductance_switch = { + SELF_INDUCTANCE_MODE_0: self.__constant_self_inductance_evaluation, + SELF_INDUCTANCE_MODE_1: self.__self_inductance_mode1, + SELF_INDUCTANCE_MODE_2: self.__self_inductance_mode2, + } + self_inductance = self_inductance_switch[mode](lmod) + # Evaluate internal inductance internal_inductance = lmod.to_numpy() / 2.0 @@ -4214,19 +4227,21 @@ def __build_electric_mass_matrix(self): """ if ( - self.operations["INDUCTANCE_MODE"] != 0 - and self.operations["INDUCTANCE_MODE"] != 1 + self.operations["INDUCTANCE_MODE"] != CONSTANT_INDUCTANCE + and self.operations["INDUCTANCE_MODE"] != ANALYTICAL_INDUCTANCE + and self.operations["INDUCTANCE_MODE"] != APPROXIMATE_INDUCTANCE ): raise ValueError( - f"{self.identifier = }\nArgument self.operations['INDUCTANCE_MODE'] should be equal to {APPROXIMATE_INDUCTANCE = } or {ANALYTICAL_INDUCTANCE = }. Current value ({self.operations['INDUCTANCE_MODE'] = }) is not allowed. Please check {self.workbook_sheet_name[2]} in file {self.workbook_name}.\n" + f"{self.identifier = }\nArgument self.operations['INDUCTANCE_MODE'] should be equal to {CONSTANT_INDUCTANCE = } or {APPROXIMATE_INDUCTANCE = } or {ANALYTICAL_INDUCTANCE = }. Current value ({self.operations['INDUCTANCE_MODE'] = }) is not allowed. Please check {self.workbook_sheet_name[2]} in file {self.workbook_name}.\n" ) inductance_switch = { + CONSTANT_INDUCTANCE: self.__constant_inductance, ANALYTICAL_INDUCTANCE: self.__inductance_analytical_calculation, APPROXIMATE_INDUCTANCE: self.__inductance_approximate_calculation, } - inductance_switch[self.operations["INDUCTANCE_MODE"]]() + inductance_switch[self.operations["INDUCTANCE_MODE"]](self.operations["SELF_INDUCTANCE_MODE"]) self.electric_mass_matrix[ : self.total_elements_current_carriers, From a68b51519f9d0ddaa0c906c9375c7b819c8b4fb7 Mon Sep 17 00:00:00 2001 From: Mortara Alessandro Date: Tue, 13 Feb 2024 14:44:46 +0100 Subject: [PATCH 12/26] chore: restored inductance evaluation methods In the functions: *__mutual_inductance *__self_inductance_mode1 *__self_inductance_mode2 The previous methods for the inductances evaluation are restored. This process was already done in a previous commit, but the commit was reverted by mistake. Class: Conductor modified: conductor.py --- source_code/conductor.py | 229 +++++++++++++++++++-------------------- 1 file changed, 113 insertions(+), 116 deletions(-) diff --git a/source_code/conductor.py b/source_code/conductor.py index 440ae31..0f04c48 100644 --- a/source_code/conductor.py +++ b/source_code/conductor.py @@ -3835,92 +3835,91 @@ def __mutual_inductance( """ jj = np.r_[ii + 1 : self.total_elements_current_carriers] - # ll = lmod[jj] - # mm = lmod[ii] - # len_jj = self.total_elements_current_carriers - (ii + 1) - # rr = dict( - # end_end=np.zeros(len_jj), - # end_start=np.zeros(len_jj), - # start_start=np.zeros(len_jj), - # start_end=np.zeros(len_jj), - # ) + ll = lmod[jj] + mm = lmod[ii] + len_jj = self.total_elements_current_carriers - (ii + 1) + rr = dict( + end_end=np.zeros(len_jj), + end_start=np.zeros(len_jj), + start_start=np.zeros(len_jj), + start_end=np.zeros(len_jj), + ) - # for key in rr.keys(): - # rr[key] = self.__vertex_to_vertex_distance(key, ii, jj) - # # End for key + for key in rr.keys(): + rr[key] = self.__vertex_to_vertex_distance(key, ii, jj) + # End for key - # # Additional parameters - # alpha2 = ( - # rr["start_end"] ** 2 - # - rr["start_start"] ** 2 - # + rr["end_start"] ** 2 - # - rr["end_end"] ** 2 - # ) + # Additional parameters + alpha2 = ( + rr["start_end"] ** 2 + - rr["start_start"] ** 2 + + rr["end_start"] ** 2 + - rr["end_end"] ** 2 + ) - # cos_eps = np.minimum(np.maximum(alpha2 / (2 * ll * mm), -1.0), 1.0) - # sin_eps = np.sin(np.arccos(cos_eps)) - - # dd = 4 * ll ** 2 * mm ** 2 - alpha2 ** 2 - # mu = ( - # ll - # * ( - # 2 * mm ** 2 * (rr["end_start"] ** 2 - rr["start_start"] ** 2 - ll ** 2) - # + alpha2 * (rr["start_end"] ** 2 - rr["start_start"] ** 2 - mm ** 2) - # ) - # / dd - # ) - # nu = ( - # mm - # * ( - # 2 * ll ** 2 * (rr["start_end"] ** 2 - rr["start_start"] ** 2 - mm ** 2) - # + alpha2 * (rr["end_start"] ** 2 - rr["start_start"] ** 2 - ll ** 2) - # ) - # / dd - # ) - # d2 = rr["start_start"] ** 2 - mu ** 2 - nu ** 2 + 2 * mu * nu * cos_eps - - # # avoid rounding for segments in a plane - # d2[d2 < abstol ** 2] = 0 - # d0 = np.sqrt(d2) - - # # solid angles - # omega = ( - # np.arctan( - # (d2 * cos_eps + (mu + ll) * (nu + mm) * sin_eps ** 2) - # / (d0 * rr["end_end"] * sin_eps) - # ) - # - np.arctan( - # (d2 * cos_eps + (mu + ll) * nu * sin_eps ** 2) - # / (d0 * rr["end_start"] * sin_eps) - # ) - # + np.arctan( - # (d2 * cos_eps + mu * nu * sin_eps ** 2) - # / (d0 * rr["start_start"] * sin_eps) - # ) - # - np.arctan( - # (d2 * cos_eps + mu * (nu + mm) * sin_eps ** 2) - # / (d0 * rr["start_end"] * sin_eps) - # ) - # ) - # omega[d0 == 0.0] = 0.0 + cos_eps = np.minimum(np.maximum(alpha2 / (2 * ll * mm), -1.0), 1.0) + sin_eps = np.sin(np.arccos(cos_eps)) - # # contribution - # pp = np.zeros((len_jj, 5), dtype=float) - # pp[:, 0] = (ll + mu) * np.arctanh(mm / (rr["end_end"] + rr["end_start"])) - # pp[:, 1] = -nu * np.arctanh(ll / (rr["end_start"] + rr["start_start"])) - # pp[:, 2] = (mm + nu) * np.arctanh(ll / (rr["end_end"] + rr["start_end"])) - # pp[:, 3] = -mu * np.arctanh(mm / (rr["start_start"] + rr["start_end"])) - # pp[:, 4] = d0 * omega / sin_eps + dd = 4 * ll ** 2 * mm ** 2 - alpha2 ** 2 + mu = ( + ll + * ( + 2 * mm ** 2 * (rr["end_start"] ** 2 - rr["start_start"] ** 2 - ll ** 2) + + alpha2 * (rr["start_end"] ** 2 - rr["start_start"] ** 2 - mm ** 2) + ) + / dd + ) + nu = ( + mm + * ( + 2 * ll ** 2 * (rr["start_end"] ** 2 - rr["start_start"] ** 2 - mm ** 2) + + alpha2 * (rr["end_start"] ** 2 - rr["start_start"] ** 2 - ll ** 2) + ) + / dd + ) + d2 = rr["start_start"] ** 2 - mu ** 2 - nu ** 2 + 2 * mu * nu * cos_eps - # # filter odd cases (e.g. consecutive segments) - # pp[np.isnan(pp)] = 0.0 - # pp[np.isinf(pp)] = 0.0 + # avoid rounding for segments in a plane + d2[d2 < abstol ** 2] = 0 + d0 = np.sqrt(d2) + + # solid angles + omega = ( + np.arctan( + (d2 * cos_eps + (mu + ll) * (nu + mm) * sin_eps ** 2) + / (d0 * rr["end_end"] * sin_eps) + ) + - np.arctan( + (d2 * cos_eps + (mu + ll) * nu * sin_eps ** 2) + / (d0 * rr["end_start"] * sin_eps) + ) + + np.arctan( + (d2 * cos_eps + mu * nu * sin_eps ** 2) + / (d0 * rr["start_start"] * sin_eps) + ) + - np.arctan( + (d2 * cos_eps + mu * (nu + mm) * sin_eps ** 2) + / (d0 * rr["start_end"] * sin_eps) + ) + ) + omega[d0 == 0.0] = 0.0 + + # contribution + pp = np.zeros((len_jj, 5), dtype=float) + pp[:, 0] = (ll + mu) * np.arctanh(mm / (rr["end_end"] + rr["end_start"])) + pp[:, 1] = -nu * np.arctanh(ll / (rr["end_start"] + rr["start_start"])) + pp[:, 2] = (mm + nu) * np.arctanh(ll / (rr["end_end"] + rr["start_end"])) + pp[:, 3] = -mu * np.arctanh(mm / (rr["start_start"] + rr["start_end"])) + pp[:, 4] = d0 * omega / sin_eps + + # filter odd cases (e.g. consecutive segments) + pp[np.isnan(pp)] = 0.0 + pp[np.isinf(pp)] = 0.0 # Mutual inductances matrix[ii, jj] = ( - 5.0e-8 # Mutual inductance imposed by Zappatore input file - # 2 * cos_eps * (pp[:, 0] + pp[:, 1] + pp[:, 2] + pp[:, 3]) - # - cos_eps * pp[:, 4] + 2 * cos_eps * (pp[:, 0] + pp[:, 1] + pp[:, 2] + pp[:, 3]) + - cos_eps * pp[:, 4] ) return matrix @@ -3995,24 +3994,23 @@ def __self_inductance_mode1(self, lmod: np.ndarray) -> np.ndarray: for ii, obj in enumerate(self.inventory["StrandComponent"].collection): self_inductance[ii :: self.inventory["StrandComponent"].number] = ( - 1.0e-7 - # 2 - # * lmod[ii :: self.inventory["StrandComponent"].number] - # * ( - # np.arcsinh( - # lmod[ii :: self.inventory["StrandComponent"].number] - # / obj.radius - # ) - # - np.sqrt( - # 1.0 - # + ( - # obj.radius - # / lmod[ii :: self.inventory["StrandComponent"].number] - # ) - # ** 2 - # ) - # + obj.radius / lmod[ii :: self.inventory["StrandComponent"].number] - # ) + 2 + * lmod[ii :: self.inventory["StrandComponent"].number] + * ( + np.arcsinh( + lmod[ii :: self.inventory["StrandComponent"].number] + / obj.radius + ) + - np.sqrt( + 1.0 + + ( + obj.radius + / lmod[ii :: self.inventory["StrandComponent"].number] + ) + ** 2 + ) + + obj.radius / lmod[ii :: self.inventory["StrandComponent"].number] + ) ) return self_inductance @@ -4029,26 +4027,25 @@ def __self_inductance_mode2(self, lmod: np.ndarray) -> np.ndarray: for ii, obj in enumerate(self.inventory["StrandComponent"].collection): self_inductance[ii :: self.inventory["StrandComponent"].number] = ( - 1.0e-7 -# 2 * ( -# lmod[ii :: self.inventory["StrandComponent"].number] -# * np.log( -# ( -# lmod[ii :: self.inventory["StrandComponent"].number] -# + np.sqrt( -# lmod[ii :: self.inventory["StrandComponent"].number] ** 2 -# + obj.radius ** 2 -# ) -# ) -# / obj.radius -# ) -# - np.sqrt( -# lmod[ii :: self.inventory["StrandComponent"].number] ** 2 -# + obj.radius ** 2 -# ) -# + lmod[ii :: self.inventory["StrandComponent"].number] / 4 -# + obj.radius - ) + 2 * ( + lmod[ii :: self.inventory["StrandComponent"].number] + * np.log( + ( + lmod[ii :: self.inventory["StrandComponent"].number] + + np.sqrt( + lmod[ii :: self.inventory["StrandComponent"].number] ** 2 + + obj.radius ** 2 + ) + ) + / obj.radius + ) + - np.sqrt( + lmod[ii :: self.inventory["StrandComponent"].number] ** 2 + + obj.radius ** 2 + ) + + lmod[ii :: self.inventory["StrandComponent"].number] / 4 + + obj.radius + )) return self_inductance From 59d7c9e14a22f2be69da021fed721f8120117c9d Mon Sep 17 00:00:00 2001 From: Mortara Alessandro Date: Tue, 13 Feb 2024 14:58:57 +0100 Subject: [PATCH 13/26] feat!: new variables in file conductor_definitions.xlsx The new variables are: *MUTUAL_INDUCTANCE -> Constant value for the mutual inductance used if flag INDUCTANCE_MODE is set to 0. *SELF_INDUCTANCE -> Constant value used if SELF_INDUCTANCE_MODE is set to 0. New flag values: *INDUCTANCE_MODE = 0 -> constant value from variable MUTUAL_INDUCTANCE *SELF_INDUCTANCE_MODE = 0 -> constant value from variable SELF_INDUCTANCE BREAKING CHANGES *New variables MUTUAL_INDUCTANCE, SELF_INDUCTANCE *SELF_INDUCTANCE_MODE flag is used whatever is the flag INDUCTANCE_MODE *Default value for flag SELF_INDUCTANCE_MODE is changed to 1 modified: input_files/input_file_template/template_conductor_definition.xlsx --- .../template_conductor_definition.xlsx | Bin 19503 -> 19403 bytes 1 file changed, 0 insertions(+), 0 deletions(-) diff --git a/source_code/input_files/input_file_template/template_conductor_definition.xlsx b/source_code/input_files/input_file_template/template_conductor_definition.xlsx index fd49f5125838c7c4d62f5a5d3bc5dc678e2d0f73..5775bc47de57983ad9fac858a072e5a87cadf341 100644 GIT binary patch delta 10459 zcmZ9yWmsHG6E%vvCWAxp;4rwmySux)1{+|4OMpRw27*IycXxMp4G`S*%Q;uhz3<*X zc31a$s%7_6)vLNHzeA6vLRTRpLUg4uYGF}LVL;D8bd)BrP*5L`prA0JprCvlSiGIx z?9H5=?U{WX9e!yVDz5?o0YB>Au^t7fCbU;0f))yGWQfQX99fd-7hNMmKUz14M^C*y zSE#F@tW$B})_XY|e){2e-h6lYyYRyViF*1dpZLR&G)7qoqfY=Nx3ZSqJA&9z#t@6# zJ_VHPS2+E9_}P2VzQie2ZG#U==6^27V*rigjT#G!Yhn*X&!?)jRK!b ztbKCMHN3`BMAzG*U2#K$hoLE&zKB%edQe1;VZ{h%bwlEOi5nmKAXF;>JIXFT`LM4aXnoPt@MU0fm2xSo zQ>5-BtEDV$cwh3(ldZLF!wyc#-%nw2dg#9 z+J$G|U86mUBvy>aQY>YY(&)sT$(6RAzT41NGJjj(EuP4(U)0Zffq`tJ0N}|KtpJJ; zCNe5elD=~W5I^jk=8v$?f)hBHV>~LT)K;q;Prs;iotb73zbb^aHVabtdLIxbQd0o= zNpS(ad-$sKxnpAlrM#I-Y#tdyfRxgRJvp*P=HhyG{RmcSV1_R*x%5kj_uAIU&7TaP zyqlS%OntboD^yK&4UtH!zxEOdl{wqOk_~J@vP=q|q6&Gu<130a6uIZdA7mwC`Ni@l zO^XHaVkU9)WsO!betwZ;T5^_~{+|7LnU_cQ)@?rKdm=g+?#~Zii(;&aj4dp71;2m7 zbT)RwTiTRuV{8-Lv;BGgNKN-$OKN`smhHELkmq}p0rm*OOv$^oX$b1iN4+MvPb%9d z3LuovunYoS1C?J5F@Si?B;^t*TgGp1grcNP1Ei0^_Hh-PJhAN11YSnuB4v z{C9rRG7-lk`|QY(k9fH|;m?uTeA({P9RO*-BiSSX&VZum(~hgo5?1o~h9-Sl1N4sp z?`EClt5D?5+#;LK_n9;q4RhK{k40^48n2%-eRHMYh4NuqmK4?nw2N}jzC zNWqu`b`d(+8$$kR!ka>Z0&hc-{o}3Q&IehCXMF;`s2I90x$0^Vl+h3nJ;4^z?$p)%lissn*J-2yQY@hg8B0p0kgrRld=>= za7XAKD?>nIpbDWe6F9;EbQKd#mdYPLDj@Gn7qE4KTD1YaqZp5M^4ku!8X?J1u&fd& zE5!|p=gdezLMjYcMM+^?OJD(-Ic%y8*t~GIwOX4g0(}~leq$y^=q5x!;TMiR4sxl` z74NBJny@_UBc=0c_U2;?MXogw`yAlGZZq7&)D(?ZI9%ny2C8s@9=-8bnGJN_vRkhW z7?)W}>M2f4h^h%<-_W#A1068+0*c&^Z(F_I8Y^!uVYfd)EBEME6NZ%?@(!m|33P8V zHI6bd+rr&kfEX~oWCicz2t-aw%l^r1XnBo0gd#p3XE%0>4{oF~KcrU}bwoJRX2WZL`#xc5-jc7z+v_Wo>~c{e8jj?84u{{}cl^ zgxD&xjaV?HCW%gqbfP@Fht#E|)TpKeOu8+SW~5Q&@RS{f{Fjh2x1-ab=(?({iP1zl zb_To#)O3&titM6c7Bj4bJ>tJS% z74q85S~BrkS912%a3dM|W1JPWs6sqHk*eb6z6Q7KkVlvX9^vGAJmGbCGZKUB^DzR&B$5n36&LN+rf0B zh`FG2@GdR6W$X|t`^W8^dcjmI#a*aDRkIe)W_rE(14B4xhXg)i(Qc&ZRt>ZSSSuVX z$}p5XWRGeQ_|n*)*DXS_4UGGOc^c)CESxwoEn>u7C=uAAq0=}!983`p9hYLZ`)HMw zRA?2jtm1e>*JCrRz7mfWYp1=dwE zOq|}5p1Fp4O?=8ioQ4>;eABteDpMFB-pti0Jw7N&hgbcb983CxINV?_&_rlL}>LWa9^~(_Q@P( z%EE5@YmvQG3%dbB^i3!H==yN7=C{|8D-ip1zZViN)(~speidLZ>%kvShOq_xM6_|< z;>IR~`)o`tPq7I#9{E|kN-T%9^O?I~tdP=%$M01j3W#DD5-Un9t z>WWo?E%9ph8HQ|=tP7FHPd*R$B|cBfrLs%rOOn5O%!thh=H0fh+33+1d{!8mmWw82 z?;8r`RDHBq0yew9IpRq90=G z(mY!I6!tf~^EYIm8(E~K=FL80X_^yO0=<~NNJI^~MYGNak+aY|rOd1kfq)}6(d-&` z!x)DbiV9I^st-C2F>cx>bqC;@xhh3>8w=WD&u$=ni%UF@5=9%?v`{ZE4i&r-cIeBf zV+{jC1LLt+YVM9}=ZPfMDs(@+?8m1bSx=^Dsm6>!ks}^B{~R9kJh;v zl|-#7D)(Vl2uMV6*znQ?LP7(~D>%Mp*)R^$UXYuj3#ctS%MZKRGe664`TV7KpUwB_ zrJgMWr{4`H)|Yz-zqU^yhYQl69SDMKp9LpT5xtI{w!5z1XNz_{VMwe>>tF^yz77zE zrj!%SvT{!~!U-qSoO4EL-ShNg_ha?VaLRdg?tY(eO|R*8IV2SW)jT6rJ3X0PKkbs# z(grOl9vtJ4e~Wp;jZl!u->YN922Ae-!SBN)PLFhcaXJ;GvYJcF>Zq(SxJtW_&MkM1 zxubg(%MR=vwvs0|#5?XqK?Sy6U#cLI43+KP^R^KR+|;aVh9MH6f=LW0v;x20yn0N$1xWW zYqN|fv6Axr2nZg02{zF!Bs+2`Oi}An^3WLSPbE0ZNee=9o3x*BJdkZmRsak4b=>d!}ae+1YRnw_AXL30J4qa>04`6r)ZW)fURI?%-P&> z4iS^Rm@TYPP5QO4CaaG%@d;r?2_gpO6i8@oVr&f9B~l!1NV7UPL57Od18aBWtsaI@xaulBkjo<{n=#w$TE&cU+h|rGaMOlq zHSWi+75@f3*apq`Z_mEUu8s)op9`R}uY%H$$NIf~1Q%af=FfT|mik9tGAkM#flCAM zZ%i1pMga)uX*nel%S0;ozmIAgZ8+kV+P&@$$LMEUJqjWRX${o zAEn?SOz5|DmCQxaGNsE-NpffQ7JC(ka7O(U;O{~4n8tF@-qf22cXoxdRzuBUM*QG0 z{m~l8q^{9T2(LnF`+>?4ySTx{xG}V*+Fai3>mNhI*GE^W@}QPaGEA7T$?C6M}mbQ@>X2miq9{hMPu3Arh#A%76OO^E> zQsKhT+m^#)%we9HAA?DbjzMk;?}-0ey<%ixi92{Gs63p12?p!G1ViWlLoh;eVbw;H zjCXt*5C^@c(Cqc7l(Q!o8`!5x-s$ zhb~a=+ji|#E9yC}{;ZzV5^W<`^&TB*}7Lg`4f)z3=G6#Nt|n-0QuBM65OAO zPrWL%WIGGGKe%s$tT$sFDsKSrLe_L3{b*bu|HxtKi~Hf)dxu<@6%XF%qu$6>n<^XoKSwH%#tvR*Md-_?!!uaU88oA7Lued9#yOzN2pQ|AS z-ATq74InPlKW>#gw!DX9sEN*WtIuMqO;xOejfq+iyCMD;YuBFo?Vwt(qT&W#l9>2> za8cn&&}CWs@i}L^w1ZUILd8l|*si61e zcy%p~naU&{U{K=Yfj*}7(Qb>??-@>TP>GLug5~_{bc=DRkB_g$cfRTfJM~R!Jtt*0 zl}nb{D1V+pCf#4(85`LS2IKYsIcXXo1fV(PndocpqO3B}FLa7SSWw_V18M=)+k z98q(7TDDMI`f-PvHS;?t*(!V5fL)Yi$-bIaW^{uQfX*259pvUKBP1Mt(o}EoQvGt< zJAX^`TpQr)%O*%(gq0S+V?|nV=v})Xo-z3)Dg_2C@3z}Trc+^~yK^7Is*SW3T2o>w zV#OIRnr2x+7K)rpxW#_n^kX^cW~*ma2Qy^U*G-G}BMV@sIGZq{L|RlvG~SL5nXe2F z4%NV26p0yi1f-ix7Rp0kj9`YM7XLG+RF#FBUJ}XsY5-fzMHh>gTTMwcj=nVbdlVqQ zxY%xPw3V5zXIIi(5{a4NxDH(s#Z5yMy-y>bOky724Z4!>_PRno%EZ~e0;yXnhA)f) zp11U`YI?PJvzSt;SJ)F`hrr6^ZEQ!I($h;6ca5Y|kf%c`u#i5M4xx15G^b)wlex4e zNMZ-{g|sE=Xbd_Kw(n2W-YiSMojEeAs^Gy;ansQzX=JSv1Pw! z{ScWUNZk$02IeSc=P0CVw)0ZCDzifu&`n3}#127?<^42*a4UtbJh;S#$N3#BW~>SB zf)A(jFiFwm#5Y@rBy{&@)*6iq?gryAANXovmok5xi&Xi9A?!vwFs~WW^4d#;hD7@7 z@rfVR8#^1X7Pm=K#Td*}Kgpv6RMBE+LdCALDbQmpOLA?C4q`~o zGD#Jsr6nHOvi~5!R^2BaPZi=9)weoz|7laCQ46UN=mUnHN7L3Oz9Mn*vS6Tsw0JCo zR(t*wQaN>ZRT9L&+`$-goHOCH_L&I;iAS|n5w9Zi477`gV2Zbq>}oSKbMYqU-vB>8KJXb5=XW{NMU zPnKrQ@KA#jY@*AVKPAfTygj{QSn|eMD9pTSnIGwb!$nPxzzXL&5M!OT8&FdP*AImVJ+ zf!}xp&G?{MDYU=Lh&a6XfmkX zJ}V2~)5~mf1|sg~lpEBh*5*#T7U``KZ2jLJQWTs;2SXbZh9mpgsq zbaR??vI6n#NIYOLJQ9R$-~&p18~X9ItNMC<_-(^RAcHG)vYlk~ZG5m=*zWXBNV$Eu zm>X$=CBWH5m8Rj*fB*LG#~<{HO5a1jNijUZ&hEF2LL+z6lW6){H*!z$O|9OvJ(`p(z@c=KHg1ECUr<26$ z=@UfQ?$_DID6NM3CkRcbIc+l6M4ma{Jy~}%Qs)LLNdjWgl;QEWs-^I}z=8xA{>b54 z=G(CIsUh*ldByW%C%3bCZK-xsb0)@^arF_nOc3bLm}$18q>cM}Au2Ck zQ(qNahUL}$E~dHZ>&h*TJG}LwABK#OPMYgNu8M^R!*6m{Ohf?kgTKug8MdzbV&Tae zzGvA7OKp8VfQ1PimOH@oVJU5l(NrN}?w1o%vyJ(^8bX$RuXjWujB`B9#x?DBnJ33* zm^5a9lv=QHJ{G7HgIOmgk&Peq+SI7gHc9<6_KJ!Y~Y5KW3%G=N-kLpDhX(SUiKj}pxO7Zs8azbGKAY= zapp${xh7~*{fS$xsFie}$A*vu*Jv3^N=2cWHuY7-V2C*vpa}ixAFPD^%RFkTdN|$*lrTBn*#w3s!i)m3sd# zSKeXCD0-0Se0*?fdUMEPyBH3BU_V%<=|Eq#580VSo|bX{#A`k6hwqAd(ehGL;1Qoz zFz$4P<=wedm(6xaDNIwF2(YW-0hRx|*kr#sJqQk(h*Ai!$Kcc{LzJQ#YLBy5-Xksa zmGxON!zpNB=vc|%kp5$C5M~lY!G3G<76@`;23#Qf7MK0fbKZk%5at)pCm^*U zch6jmCG6)8M&J?Ul%yb>bG{ZxpeQr;ZJtO*)eu1g>I3jA&}TSz>RM{#FlFFkJi2yE znB)6640|Fi#V>zoSnNI}YbTQDbR0KuNrL(n93`t%FOHZCWO5^rl0$u1szbR`jPBrV zNkBV)Z7{)65B#`b?Mxv7rDryDP>zN+t$-#GftJMCo?V2xTSxe?;@TT9x>|kAcyD@S zZXZF(!BQ5R7@l6dWbq(olXleGQ988=u7}yk*DP1_#_jk>J3#_!pG7d_$k?J;)NAbv z1?<**|J9q|b>*P?4Y?3uQ^xXM&EkE?C`hz&dRzm#>WpIx$AwBs@gY}7bClQob7*Rk;?1`&M4xqT{~;HZvDH)MxXa>xCHJA+n)g22{^x*RPZ?uorZZn z-`U5xE>o-#nN~l`fCRD)9`ZKoy6;5DeY(H;Y^|`484?AJ1;JE1{JC?vr=O!$f^={o zfz(_Cclh|Bvfvk^172uwQJG#i|jxZUHWat zw}Wr*U(hiXQfi3aC zJvpcC2C7keK|krmXl6}5WPSY1+>R|~0^;egS?WE9+dUC}WFz4|B!aCJ&U$#=>L*zW z8zgn^Cpy*Dg|v^kM1`@1cdxIpCaThK!|X= z1Fnj2Arf3wm)oj`P}Cppm(SK6Ehm++{W3C${>JWhw0gm<3A^cQ&7ZaIaz=cNe6fTAmn28EOMREW&fA|7fV z()!_38SfTV+z6KL6qJ*WcSc1$C{f5rQV&3nOggW3jgsB;DQ_5^HPZ8%CV)oYdLAlz zG~65*;!Cu?;#*p8RjxpTtAl#!1?F!fP536#FseFjp z?s}!IK3TUynGi?g=;yjBY50)qg!LNUHNY?B>9u-!%jp$cFQJpGw=er9gZ?KYq|;bv z^_Gijf>4OE(gk*F*Q@Dj+0&63zeKLmKQ=kkOJ9kZM~nw6m2$u ztPzsCO>dduhMD}56aA>C{yFl;OMbjhiPkR_xWs~LhOE>diH+d!Rntl`hb2SaTjblJ zrIN|nXbfc2L0dTG7p=qRov{qHJ2!XPYG5caE#waFLgL53Mxcw$Y;``Q}^(uAGnby6aK(ZL!on8a!U-cPJOCq)s5Z0nlXpM2)`@7BJnTp() z;C42G!x+XSsx`hq!8}JO2Pxjp`sSVU9b_*<5Lq*z1V|e6$%2H4?d0=?T?V?*wdrN6 z4Z0U1*D>@V>7b-%weFUFtvJ~OTiLjS#7JW0DOW8~@uTnLvKYz=&bp1Ija9k_s=ui- zr}$L+%qgluV+I?T7vP3#Vbxm`2N$E_goUNO4EC&aZ!yM0EZJ0fOW(BB_A`eSJzFIQ zU39oz78QI~&uVf1Q{;Nfk~gGOVserIN~y|n;$YHwWp`>-+K2!iTh(e`hY6~7FIV=4zak>F-YZ~K}QheUcc z1!mI3F10Syp7fPalx56x5J@`6kG9}}*!^$(9o~D^DlX--g8luc@RPHWIg#?@ar0@8&v6r~QYU zOs_(v1xyZZ8yL5Yst0h-4Z2NgQ3R0P^x~fem2xL|cyBxwx|`Hg~Bi_Y%{jqKC(i$fjyFIB1?6Q+O$ z8O!$q0sWiso#MlO8J|6W=Ki_B^m8NQQT0b1Ue|I+OtkQNzAOBr#Sb8LKW`SE;HCO; zu1OLQ?3p|bPazvcG1W69=ZtdGqSwbw>y~XlY@e4ag6LvnWMs7W)3is&r=bw~guz9X z%yv*hJ*WI0y>}^+iS{fpXUephOUsoMl~2K=Rkh$(H6&}pz`rhv@QPq5W6yiwUvjp@tyK zVL`lQ31I@TAeP+!ip2k=VEnHl5cXe1I1qdp9OC~jkNy|*j{h&n=r8EMSnL0Sa6bJD znvupK{{OS@{{{#JB}(;|i+?8mU*Gi~gog=&D@*!!(-OkiutG2?aUjPsSdb1mLYPAy z2oXOOLu$$nPT6*&lEKfEx+`fCT^mydBuRoLueAoSf|0 zyd51XH0_--1o2;9lCNOFAz4N8rC($FC0XahxhhzghEga5ZdW z`W=cJ%)WQ-7MtVdcjKa))Q7eaROf4whu6<7lU18Js|A+7K_Ahk z9*NtdcoEPCzX+NiEdT=a&_~3E=<6Kgfd$>0Av^`4Y45gx$rC%uO*VX8t&GU=nzr9! z@fLwajSEx~PlV{aIaATJEcPI^VayOHuKZVa)WvWfy(BU?ZPS-!Uah_|L`tiMf(y;y zbg3lC6lI<*Z>^vud3-{7%pcg0mw~`O!&7*5Bz;ZpErps^K^U8FH2&+tt_~qLVl{%% zT`n>mql0Lg)-9&JJzjOFZ*?}_pudB}go^0tfpFkYvuwITGwic-Mr~!u1rCT?Y^z)f zmCCqdaM99sH4^+xvc*zZ2UY=28K&~`sUC3KHwuq_G<*6Vwk=|p)M7;;1NPnn(cq(R zfqE$2EsS&_}Fm zisbMu!htc+#pd$)Ia+dJY`B1~Av2k>z-6l_)+T+4zOJUfaL zGys%?7Si)~;HVdoMLu=Nc^#EWmio%=co*E7cSG`NK<_O|0d47R8cTm6DC^^l?u!O1~B8op^ z^!60|+_D6JZe@H^0;-&A#%uiw#qFeNB>FKD?|S0M#ac?nwX@bY_S{IejB7A=VyNF> z{`+~6*W6LPo2Stng6cJ6rEP5Qbj+%I7|nV*km=mAjx4cm;$UBaxAJM2vBwgTN36Yp z_;JAH(GH~Zq$~3rdR5>oIN-&ml-1$%u9QIMQ}wCd`!%=9dF5qpK1HZIA)|lUBi92# zb9fc*uLmI&(*K4wG!pdJXt!XSch z{I$3X67`)H!T9}a495ZzpBN& zT@iQmreYgF(t*Rf*kjOk<<1(J`TggLbpV%#xF6v|eU60?2wrPm5@^Hd<|36^))hD@i)I{>u>!Cg%DtRm5Gy&`D53J$a!RoiSh$Cl{pzT13f5*)NopgGFGcns7)hO{WicD4YOox|R|(`B?z`c5#MbnW$X)2MWYYuc!gU~608 zlM?pK`$bNEdS)JBzkrl0C#y<4hCHuTVxzCAW7sKdAjsXL;!5o>8oyBLCb46_zg*m$ zzpv91uwzOThE8AmMox20(p39hS_qUyUK;8m_0aRN@gA9tE>&*%8Zu|DcU0Bo{JlE$ z3_O@AU&Xw=MJTgjBV^^%oH$s5ax40i0f(XQT4_D$bhA+oEX-7DVX`uh}de&+ne{1@lP} zJJWP-MaE&VD2=bwx^pj&3z;cDt@~+>%Z9#xTOMAxJDcpM9SmZODq6k2n&VeSrqg;L zB~!myo2Q!vwMT58pEgsf@-z3i&@B3?vVq-sLOG7Zd|<}K%8pSM*5ho$MVL%K<|BBdW-2O~sS zS|l}VVAZ0Ok830=Z?P7n%ja^lvJdH#mYmes!Y*$!ve=!@k)B&SGFVRuHe8=S%)0(8 ze|nR3`wLQC{F7fQ44!AhF5M^F%5B-&wE@Djs2w%1nyre{HN!2MZRInBDCj*3 z%0pE^)SZl1ZSnZfAmBN-aeJLIW41uAKY7q##ky#kb-#uQX+-CJ@%3nCZ7P{g&7hCj zI17oEBUzI?B}LY9bm7L{wEq&HNe3u!r!)zC9)(P`1q(`Dn1(gKi0;-Lsav#@UcbYZ zqQD692jtRG>BNVG`sRw_QDF{?Norxkb~6@;)?iawQt=G11{Gk|DjNI{Cr$pm&${Ru zKp4Y~%nMuKP=MIln@>mMNqI=(qQiT5`vr$miUuPn29Vou2df>7`1fS<%-c=yp(+K0 z2*!+il5!sTi)-BUSJB>E(r-`c%V(KVN<(B=(W-|R(zF4%P{_B8cJAm+x}wxc#OZ>nD?9G>n)iW^JL~Fng!Kg59JlfcQkLvI6_ZhG+t&<= zEw{tZc@9_c)`})0M#NGF1mH;_jUHx5`vm3rzz4g<)4GprFHo!Eo;{0KXWei72Ib<2 z?Oj7Wedv1TW7yFJ6``8UZ_{~7_#LiP(RPD4n(O+$=MGLu&smk2GE0(-QfiQiNkx_n zsR%Me3hvgI#FEi)pNFIiMwS>pQ9mL5555gl^ijLP*x%&_VEiblWgx%fC?BZ8eFJUg zJUK=~k7qZJ71cn!KSscDqZWz?lnnHN-zeW*Q&z<|bS=`)alJRGb3E=FNoym|vk$No zz`YvDF?qe|e6v#T1vC*gIXXIt#(vvDQ&`8OYU64zC2P3KU%bEV-Hy+%XWA_6yM=P^+4EC^HV;jj6&!iqVIN$|}jHM+ie<0}pe z=B$z7x|J+$aa!v``HtuT3)ZZW(T(tfNbpYA;_oP>B%j9b4eV0HmJE;x>qheMA`b^T z$e%#wzYI#Ij;Znz82cGy7(!34+2EK)gm!J5QG((hYm5;ZNQdQlt#zsMYIU05e$3e!jbfDe8*+>Qf)Z z9>(l@x|Or(<}^?6Oa&3GgUZR@YS=l4X@v4J2=Dz4yN9fpuj-HHxXM@;s|4Tzvb**Li?_;$fR*7W++5(Z*u;0RS@Y2*S4gBQ`)pN*S~YP&;RK0e0{ii z0dB!~c!tw-^{1Or$|;k&nGQi+Fl53ILZjbh)foJ+bq)zK`M72F+-Bg#5fch+AqA|>oxA(N>x`$l(vz7AfNr^a2;`kGArl@-8mf@C%uspWN zQl$`mrKz2Kio#7?^B$F+a4?r&NbiHn|P)c1XNu}{|;iM; zR#2#9F<&N8zyk_Iu4JlK7v`WI!05cg;XG# zAFn>Z`l-K#V`+s%o^k6dY!Sp@1&hy_U6p$F2#9^$fK7_|Z8I*LY?8~Lu%pEsZywtL8jAR?e9n^opjLmN}UVBB*< zDB}HCnLod*y%)0mk)?o27Pf`eZchzF)j>4m3LYtueX_;QXBQi*Y-t+uk85e5vLpw` zf&H8WrkK!^h@8ACqpPZV|CW z<(tl5+iAh^TB&#PM{se0Wxrg%$<9xbH;&e_O_Y|rv3DMF^zzRr7!nE)h)lcEp2p#w z@pdupNWW!dB6_+fAL!f0-915ip~8w=mq^PqT^JP~;sGwWuHwBVziQM4o?lS=R+WE% z!2Gq*!2wj{;o$M1S%N628iUA=QSbmjJm%kK8BBx-nUaso0AmBzKNh|q4SRvdMN|gL z6vNw^?PmHHIOuh{QWwYMHD5m>?bU*I5)322N(w^b^F2ogi~Qhk#hMOlYF5z5Q!#rJ#4^}&<+Q?4IyklUGD?2lEj&2)z8cshu^vAAL; z6ZGyNWp96?r%x?D%X-qH2BE}jT-~wYKj4V-#F+S%&MRH4X5)+s*ik{IBD)ND2FW*Q zXO>phIp~lmXY=rGwpF5OLD_HHC*4;>SNg%Z=wv(wIDVVIIZtlCB^8+H5zXjOSh$Ye zcRVqr;SPP6^nb?UFcsroi#i7wm+qT0eLfh4L>krJ1R9cGo7c@lyn22_ai(L>|I#9ql7Wf zEn0R0N`gK7Fxa79nK4;`j3Muh&#FM@oanj5uY4Y)aIznEm(FL5TghLlR7XkLeji0a z`2Dubru4i`)k+*3V)BGoTTLgRH4PMn`9m5=@XPV|Tf+$NCF9u0B^;M<`_pnHY?H&z z1ohm^wavIW%mmDnp8P)zS$eCl8)U+s=h1$8gE+w|C)`4J z_+*Ej#^%CKWO4YTd1WO4#IR0a-aVCEd`xyZ8Am(V%^2n9)O^qB`Ieq5#+{^6{myBc zE5>OyX7ksQa%zD~-Pi|($H{}V4gwkhDZ)wqCyeR)$qnYodQIBIq=y|j%3@OZu1;g`du>%FhD5^(DCg6YrZJFoLP?%S33YH|I&p{v ziY1*Sv>^yhDw-Sk$TKlMJ;PG7;%#h>-$;h5KZq`0uDxJFCf?6Xh(sT-78)EaKbt8Z z8gblsyxGdB$LH+cfQ~l@$pl9Zv~8Q?A{Ti!5Qmfv$?g2Y9BA*P+nOcQ*!StJU*pq= zrp|_&+!p~iW@GK}V<^4ucWjCU~&Q)urFOa;ugsBw>hcZ)kNON_37_@&z>X zLzTU!qq}bf8%f!j$=TtATUQi(NzO5b-P&JTK;KWa z#h`^i6_vNZ*rBH~j}$CxBn-A_@`?DhV#psr@ZYTG(zvyraU8;$^v>6aij6CVa=Tpx zV=H~E`aT)1Ry0yOt$%MQ%F|6ygO!@sOP7>1&OVhEW;zedNC(+}%_BLxvbx#IGth0m zVilp30e9p%A{%z460q8w;4KAqW%coo!o5J+5H7RejIfmSY6TV9({?8(SaY5Czz$O~ zFkV_$#0v9bq>SWX$lB1@sjXgS@YFKv;|JP`6KswYxzJHSp+wZCi)qD{+h>pJdK3bS z}TOaKT_0q-`7Ac zXdIpt=XYS<#0>-6D?ktID~&jdESLw3qnl@??NhqC^4~XffG!gnSG=H4JBg!1qQVf` zil^dt5w4l~-=az42 zB9kPT@UM3fiDE_gfbxdP2imWcOEV2#tOhfG2jCGov8fuJeG9w_wx{ByT zd}ee-Txg{lv!PlPRQ*!UAqo8|S|RDkAwkqGDw{=e}0Z@AW*awH<5B5}S~d;Swr+yJ2}tU$=phb2C2f8(#O zB4>ubdXwjPSqMBxVG!CR11tjl3u%YwyQAcOl`i+bJYv!f6m;A8 zk~#P&=j^WRF_;)#1KBW@uq@-GC%}1mPq-#~ zt=ZRvCb%4%t~BJdpY8nRZazy%tv5*u(EusAoGt-*F}j5|lP-P1^c_xgHHnJ}LSyri zOHkl(+j-xg7k`K%4G;dtdwuKf`%Zaap{UW!Wmr{>(d*0TY=OO8_9ZPtOp+7u$pT%#i|!L zGQ2{r0l19)+Q7;1B0MEGxRP|%eM;$3+UJn4dEd(iUnOXkcjY z_2GYN`dY8(3#?RL1W!9LLdH)U%G+xLDQ}N)yw_UqF6C%7`mbAWhFcz*RV}JDYAy#- zfDXm>if&qv=g6fh=v&z6wuPgG@21L$6eXAG(Jq-IXL;?u4L?2f(1au}b?vjZDEUS_ zosy==Yo&~Nm_FEQ=7cc_H)BSP2f7>_bS`w&q+{>r`*E)aI8OJ_V59w?xfCNv&kLG_VPq#8AbFqc-;au_6A4)T zp(sOX%XfBLxmEa?@Hyj*8y^HWPEebsUjH(S&?G?|)DQEwL;lDB03RV#mzD&gqwloD z_vVwa$!m|UaHGf}U#~lugn@`JcCzHmk?f=1kr@>j*Rgh{S6odt(^z;dmxUqzsOb2r zD^wyrZ0{j!c<}VrR;*c#B+0tCGT*aUk`wJx#&URS*3a{PxA_ebs7{8~hhACV-d68N z&hsS&L!VCf!boP5nZ0&@SQtGd+LBp8z<3bd_)V+Q=iYK7^7N!+YR=k3UPa6`PyQd| zMPin$sF&3BEw|bR>K2FNVmHG{YcDJO0vv}&KeG!N!fuYXzZegFe(`qaI34esdI&S} zH;oUs^dFrjsIpd9$xv?ZB=Ht^;{LY!BVRaM8S!P`-@pH{6oO;mo$@jO$+_QBeZDNl zTD1|(;Yk^9CmVjb4yzTmJGv1X{eDpL7I}=_-^p2(q3PLg=lUk+7=vEB{j%LUH!6G0 zIjc2JHsBg#Q$t9(mxb5|7W7+=Vd*M`$axvI0 zD{+RJtK2{z6|K_JN>-;1!kpwTmvzD>O`V}mOTY}ZFqd_<4OWRyWM6*VBSp30n~O0h zV+n>WbHqZXi7c$)%*){MAB&PmFkrE;k+qCX7{awZN4Ub2VU@A?%`eOu^@**EVuPz$ zIBi_>oKy1WZ~>95bWM7zI&U-+dy-DY34S6z6P}&zOmK||-)xNx1h;rcQmrI>JLmJG znM{aj718bdn8P}F%+!R~x4ztkW^or9_fk^Z_K`hbWz}evkUe zmBC$tJX6f8RJhxWBKG046zUNz{dn6n$xG0V5v`iIcuU-8gLMPDhr$#3JtKe1X%gPq z?(TWPLWR!zw_~Kc5M7NB-Z{s%FKB`C=<+p&^bp|}@SEA{6@yw8XU~qkfbI4qB5N{K zNvXF9CwOiUk@$e|u-!Mf@jT6r>DzrJBuR|lKIMsvA|}IfNaKXYWetv?_TkHp9)^C2 zYVu^aoM`3s3u9zDLu6WiABWqqC}B?huA?0RHKy|cv4 z#jxz%Y?46jiCh_atd?RC$K%J7t(j0mG+BS61L$4H=}hFiT@I6$qDLnkGo1GuYxkuC069!ck--J+ja0)jQjo=^u*)i*lO z)>Dr!L%VbI$2H;^1Goy&Tyz}XYK{=6Kn)5?a5esvQN2rb8NF9P8i5-Gg5e?d9Xyhm zXS}r`-D({l=CbZvi{E>tSH}DfqpwA62=`>IhrETALL<-scJ3eq!FcuxADN)GW5<;~ z;lMjY->H15ZnPKedwayHvk&Ci*&#-Nn@T&3U>Uf_RwPu*5SGFTFF>KPLuu1Xw^{N<*`4 z^E9IpPA*&hccIQHIp_DA`~om?w#V$;V!$A3H1PB<%?tjcOheS~2&P zK#69QE&)?_?e+t!T@L~;p$4C}iI&K* zQ%Xv0#0<^tuDU2uQGvQN)QLTuYrlkcZ* zT(*MEZJ_d@ZO^Ok!+uMQtGtq3!%v@*W-#aRYyxJhst<7M*|b1eJuC0*a@<^|ecqc24mYp&%x=R@(EZPg(5xdTB<$~$GGIUnBmzN?BQ~3!ib+m`fu^oQB z)T}nuO(*rmNM-n1UhfiuTCILDbIzQ>0mcrbmQjArqGaeeM3t&z_A zDs^rWFj?`XIJ^|Px3I1g9MqQ)6-nQAEa+RvmBSdK^Xm2Q55=( z2=RjssNNZVrdk_TXIyQ1qEh}ydxTQ#r{@7id6+(ho{fU0H#k5NGsY3FgR!oT$P|Ai zmLeV%zmH{~X2x;BC>-G}lNP>eyk{vnFmBtsY*UXXH>sryt)&*U?}eO%L7z@JStJaP zdc&Zel2tgoHQ_Fpe*I7nP}egdy6D~72}^y-c}DoF3RfS64GU8rjs*t;MU}&@M-a$` zsV5Tp2$O%!2U|}rCIEvHbMG&w^2cqltj`jrhJmQc_ki&~uOYmO4}Q#`mn_Xa7O_;u zga<`BQME1wIT9y~e0W`LSg{J*G(UNFyOk*)+1ufm$6ll4(I!DFs|v5@u6CdpVf5QF zGk?IRiq452gN{EyJnHbF3_CPpth@dgU8lNJGO&OQuSi1w6r*Q8g0*3pt}8dV<-}pw z-y5=s!E4}mfnkzf)ZBtFPIpU+)j6LTu^*ne8O<%Q6b)bkb58`_U_wLIDXNi{+Jp(D zEO#A}gi|LD2))LV&2Yz0+SRSi*KTEVVmQP;1+8f-vzyQ89BgR8s%Xk|g&>>===1nl z`vh7;H+(U?J^?Y?wJ1-qlD73eOGVhWsL0lV&Ya{;VuAGF0mcwB_0Q0Yr z=U}1ytmsg^KOetIP#XpcmI^2R&sMs>qSOC-JR(9#rKzAONSK)aP5SG)hYFpRVj=yf zH3tA7{P#RS3j@kQP7K|X!6yBu?E0TfO?XgI8D`Rd62Sk-piA=SDL^Qg2mx9si%t4> z76ujoU`O~*MqLdvNe%rXOULj}JPZI}{4@J!p6?Hl1llG|{0HU4uuAWs z^m5!#BB6ityz)bT$-RR%$`Zr+3POL$zJtoj{}~EFE#;|6{$~;VQ5WZb3DF3plc$C1 zier-gn=B{>09ZJgtGGHjyRn-(JO8E0pMOjYOCbsMl;glh_|ri3mpb@>us Date: Tue, 13 Feb 2024 15:09:21 +0100 Subject: [PATCH 14/26] fix: update temporary solution input loading in Conductor Updated the instruction that manages the odd behaviour while loading input files (value set to 1 are converted in True) in method __init__. To avoid that other values are converted unexpectedly, the explicit check for value == True is added in this commit. For instance this allows to avoid the conversion of variable INDUCTANCE_MODE = 2 into INDUCTANCE_MODE = 1 (because 2 is equivalent to true but is different from it). Class: Conductor modified: conductor.py --- source_code/conductor.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source_code/conductor.py b/source_code/conductor.py index 0f04c48..bcc22ce 100644 --- a/source_code/conductor.py +++ b/source_code/conductor.py @@ -184,7 +184,7 @@ def __init__(self: Self, simulation: object, sheetConductorsList: list, ICOND: i "INDUCTANCE_MODE", "ELECTRIC_SOLVER", ] - self.operations.update({key: 1 for key in keys if self.operations[key]}) + self.operations.update({key: 1 for key in keys if self.operations[key]==True}) self.operations.update( {key: 0 for key in keys if self.operations[key] == False} ) From c613845deb496e0d6279278decdfa421573fbe09 Mon Sep 17 00:00:00 2001 From: Mortara Alessandro Date: Tue, 13 Feb 2024 16:11:46 +0100 Subject: [PATCH 15/26] fix: ValueError in method get_electric_resistance ValueError: Voltage difference along superconductor and stabilizer must be the same. Changes in method solve_current_divider. Class: StrandMixedComponent modified: strand_mixed_component.py --- source_code/strand_mixed_component.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/source_code/strand_mixed_component.py b/source_code/strand_mixed_component.py index 350908b..218d64c 100644 --- a/source_code/strand_mixed_component.py +++ b/source_code/strand_mixed_component.py @@ -641,6 +641,19 @@ def solve_current_divider( maxiter=10, disp=False, ) + + # Tolerance on newton halley increased in case I_critical is very small, + # to avoid inaccuracies on the divider that could lead to potential + # differences between sc and stab + if min(critical_current)>1e-6: + # Default tollerance in optimize.newton method + tollerance = 1.48e-8 + else: + # Value found trial and error iteration + tollerance = 1e-12 + # Other possible solution for the correct tollerance + # tollerance = min(critical_current)/1000 + # Evaluate superconducting with Halley's method sc_current = optimize.newton( self.__sc_current_residual, @@ -648,6 +661,7 @@ def solve_current_divider( args=(psi, current), fprime=self.__d_sc_current_residual, fprime2=self.__d2_sc_current_residual, + tol = tollerance, maxiter=1000, ) From e7832b9510a70369d71150372d8a7d9dc9eb6487 Mon Sep 17 00:00:00 2001 From: DanielePlacido Date: Wed, 14 Feb 2024 11:00:20 +0100 Subject: [PATCH 16/26] feat!: new variables in file conductor_definitions.xlsx The new variables are: * MUTUAL_INDUCTANCE -> Constant value for the mutual inductance used if flag INDUCTANCE_MODE is set to 0. * SELF_INDUCTANCE -> Constant value used if SELF_INDUCTANCE_MODE is set to 0. New flag values: * INDUCTANCE_MODE = 0 -> constant value from variable MUTUAL_INDUCTANCE * SELF_INDUCTANCE_MODE = 0 -> constant value from variable SELF_INDUCTANCE This commit was already available in the pull request from Alessandro Mortara but gives conflict that was not trivial to resolve since it was not possible to visualize them with the text editor or with GitHub Desktop. So I add the same changes he proposed in this commit. BREAKING CHANGES * New variables MUTUAL_INDUCTANCE, SELF_INDUCTANCE * SELF_INDUCTANCE_MODE flag is used whatever is the flag INDUCTANCE_MODE * Default value for flag SELF_INDUCTANCE_MODE is changed to 1 modified: input_files/input_file_template/template_conductor_definition.xlsx --- .../template_conductor_definition.xlsx | Bin 20637 -> 20904 bytes 1 file changed, 0 insertions(+), 0 deletions(-) diff --git a/source_code/input_files/input_file_template/template_conductor_definition.xlsx b/source_code/input_files/input_file_template/template_conductor_definition.xlsx index fbe93ab2d5ce420fe9e62cb64928df875ca3d9cd..2f9a33fbd83cfa3f7adf5b946f3d38f4a75ba309 100644 GIT binary patch delta 12269 zcmZ8{b8u%(*KIhlZQHhOV`5Hh^ThTqm^c$oY}Zq6xowpZe~7?tS}@Q@c*D z?mDNsYxU~2ckVav??vEsa4^v7FD32cAY>R~;D8%7fMpIqy{nO38g2>Wjn!?K2g@Fk z%{h++4=p|Z8c$h4i6x+W-s>g%r$ZTmN57wXc#?1$FAvWs@N-7^y>z%m*XhRbKqHeL zT3+L6iIa$uwDN3A@qktWg3)}>n{vt9;pyWJvyMU63?I-8R*(|lSFV}cRzs}FKsH7a9Rv3oP!x8n1C%0qTT?rPsEx2hl4k-bt zq{YE9{}DF?Yuexs7$EUs7MCH@%^%qX1dHde?&424Hhrl0^Y1De2T#xYlO89nVT&!B zQbezXqp>gie_vEiLNgq@F&HAO+!Hb1%VD)>nUqKHX>eX`7~o#%2`?-hh$-z-FBn9sHU5#)ube>$U^)uT z_2N_0Qeqz?w28;cg&euY{aCP-oNryUc%}?~=N!Z;g6_t}`C6Bw5<1n9{~;HhviNdy zeaUuTmckPFNDf1JYD!f>Y(V}WdN(IP!aIK=iR47 zssAS-t#0>?`UPI*DQag|RE+x%K)vI%JEY2586`V+?w9px3?(!#w9xT5Y{lj0UCrPs zUV?8Z>UJS*8~}zlF9$msujklQ`0d3Lc7XSD$MD$&{B26ap&w6=U-szyo65B3MI?1q zRo_9zuQrS8VqU&y6ffq9w{h75FoH5yri?)ZeMi)cy1tV;x4HQF##rg+hZkG zVwB)hGRA}1&PjOiRk7(Wl%h*k0%Y=P?SO3RkI`H+Tj_@G(&h{W#wAIUTX+|UQ%FX# z?kPT7*T@;_57`KO)p!f6UrG1E}rRRq1_XIG@V6^ zjQ}{o#@{S&R-4*0Zj5rU3)1SB9p_a@A=g&qY@$${W1I0e}yK=x2#HHHW24%?s5yQryv{M+Cz*m)w%o<6p zg^c2loZYX_G7SKtU+?i+6}qk|2^O1b#vOTXpdtc({i2gFznqZTzouJfsOManj*~YADCqoi!+Z3JYW;9KkL@bap z)9j|V6+65TF0kiZ&Hap*RRB-9YnpoYhe}ElXfnT3yt4u#f|wV9Ob_L*sD2+S?_e`u ziJWJ~*TX9mk4WZPHoMUB{33@gg}pt0*(0U42CoQB)y6b;9ldwdH2rjUbfpDyl^A=) z$Z5ascPL@1uM>g;_s(RcKOXQIM*)HuxstISQ+tdN51H8%T9QO{ zlvc+x7NQ=%UnW~h#RxZo0Dr#*NdZ+;{&|0H1B`z zR`?U)sI;F6Ni-gI#weB(B+txvy2FIa+ELYiXy^v@*^Wy?GJ=U41D%~sF zELJ|g-DUw7IdQzM=|<>oCfMW$whSa$@F%>9&*Gx;Fc8?QJA&0Kv*wf}p5A997RZyr zQp=NS`QWB^^Q5aVsCwujSb&^$0nRB8+I}6QrUf(^X`DDd6au(9wt8q6x$G(G#?ZaM zkm={B?LsD)sCn5?6AIA@Y1d4IDAttGuDpV`K_Jnpc_6PotGf1bonHpF)C?Rmu4~%zYwP7%av_l z(m)G@3bg4|p-GwGKfpn!=)3UBcK9RiFkoQFME{{z(5?bB@c*)_ zWFSp6DCv`=o=c^+O;$z2Uv?#|!ynCbfOlW_+NlSD=|9eCL)a7M@BGQvr{;lC^HEW> zVLl@?6_F&aENArY9-nWmZ4E{ASx@R9O{@Bxkf2zUeU?ZR%aqi6LsD$I*A}?>eT==(G<*4&bxE_v3(+>+7MInYyRkx;v=G z6ZHck7S}t592tZh${W%==~yxL+{jS16egHSckX;9bV1|d8z*-mI$q6H=|ej0U)C|B zleE(P_d{HBUNaBRW>EK@ z#b35fe~4ZHxrg-|A=VVxvsIjK+TbSXq=~A7`Kv;@KdqVKllOz2T1KX$f*aMxt+WNw zl;hHz{_sI)dNo_RE2>ri9bi^Zz?W(ahC;7xd}T33D_H&PG~o5E ztbso9kvDJa0QJwokWcqetxm+ayspk?J>dt*B8r7-S%%=!*+seH))QQdzgj9Ww`{(P z(u|Z`4TWYWz}Sr{yjtLzSKd#v9z(%|HseITA(qJ;Es_Hao6)dH^oLKp3Qw6(x!+6@XA7}&?hKaZTRrYak(s1YCXINze>wDf7{n5yC> zvlPu1MAuE$c{DrRNyWq*e#8SJ0-m{>el-|RQ0gT4zqfN|o!HklafZ3GosPyW-d?_U_?*dCS5AvZdSHG4sS_#b`pMHp9`XhWLmoV z5Y56RbRK%P+I}{+ZG;%)ft|u)gg&L|d~#(h@)CYu3%o?^C=}R5k@ZZmvB{SQpj%}v zEnIKidIG=h2w~^`cr6cdZ-GxEzL40eB5vUvsL~3{j5i_q^70IU0&0d&TT{)L-6-t! z4)jf#ZVB$lR%JVB*#HtJwXXO-Jl^h4_iuUIkUF*~)?5ugR|Eu95cADpuv@&y(|O~^ z-D`}cZ&4MoX^K+((oG!$W1oTpT*Q+VDa2a}_zirEG-RW{X!>2Et6hk#tLE)L>|-hp z4^yC)q2==`v%27u?VF96)!;E*2DM;RTo`?!7j($dvqiHG>ng=Bx5Je9PAU_j`csiY z@+WF!<`;LB9kI$J_Jp=2-czmw>*lx^q6k&h_$T-YIBnH<)oM2P$}e5eK-y@$vq$Sqg??iS&1Q>o@}rqqGJowL?|r-^MMrXj4R0#nPwvWC$F}u)I5& z(orL9+Ugt=u;(d|G8%~6>|XNtOQz6B8W}KfU-?kAB?eK0$>2GYT4VMdjN~@bpS#tv zS4;BcC!-eobVOkW=&f-p0Mw18Pq#l%r_?+DR6K+$BblC;L?kcH8aaAP;fEobc2OWn zq(P)^hzggOEbGz#NG~>8M4c^m9c3rz=|2$lUcrCVXoVyr+}iwfO9$P!s7~i4H%U+Go5iuzhZ4bYi-H&ss_2)+J#6oBOBRYOA?GQtPk%F+u?dN zTtE{S4c7%v?EMgM&XTB8t3GsNgiJ@UzX_F``qwM0`kVcfM@&Qbm4^t1A9aLdvcKgppAs}zqD}0a8Jaz zLqkmii)Q*m+*+riqGHtUfOm*IL1Tk)uHtIj1>#(Ovrk@bL^;Rl{pK}Kh|nT5uZ>&0 znGg9DO}fjn0`r1;RHUc23zNd!Mt?!NS#VUw@p|y}B74o~q`I816&AL7+{M8m1IbUT%FOI(El|sHGKhQ(`7Ct}h_U>GU1? ztrR<#Z;`SSnx=gFl0T7d4cv?g-!xGOlVZGKRH>5FkVo~AEQ(G721y6j%UyY~w|9aP zdU4|W#Iq=M=&!TeG7EXy&ZL^NUO0>)$=@`4_)l^0#x!(*EOXU09(Ff6cTeA@=ZQVd z1mTAvQ_ABcDSa0?^N5jwB92i9yoM7yU16Q<0(n;P$F;vU{Bd0+yp7S}lVK8;Sk$=x zGQn}#uHSF7xT;g(n+ay`#kOH^X&^Oa2iIdup?^8!#L7xe_oop*PHoyc|0czE0x`2e zf|;Cr)bMLy2<9h@Y{Jx)Q7l0ys;kjU!NvqSqTF-sKnO^22(Ug>)sEkrtn+P4#~a_VaCw=cXgUexbxu+3f{0aj~T zt7}+Y&;7WM(^EY@PM$nt`3+T3DdcdW2Ya=iU5n8*hqL3tk;W7M|H-pz;qAzAlWhHE zbRT4iI+)T*gBg;(xh)Ov=wx0t3kckD@HwVLI+R|33?qDfG%Ciaez^A14IwUOvnw7> zW7f`EbIdE17b;Y&eH)v)tJ?^kcP;K<|8Z*Bx&4TJQ#-u!=dG%;)_w3ePGAVEWae4Z zvCAwv`1=wDZ)V%IRpW0POfaEe{@43-gzK^h5qpDb;!&YlOcI0&VUk@7b$}6E4(Sxs zU>}eDd$3jqDU|+zLN@ppiSiFclv9h|;&_48i zHidnk?&X8H5z|5m$y8JPXnImwg&_twNi`0$2-Ct}Q>J-(TB_F<6$_+_%Q6;xBegy% z4R!FL9+tQp6^Aa#v*rn8Q$SC;RyD%wT}OCm$OxP_s6y&`U@+Gl7ejl;Cl318jM z2`<^xSiT89$ZvTF#L6E(hfN=;HU!hHv^MS*B%?ZX{(SEmr9spvyCh>!w1n+1{grmQ z6i~%JxcZjzTyPXi!t(AtybskyadP=|!D0Bx?A=-i1KD=TD7=3^Kh_+MDT$ah6~5bI zf><*n!o)3P+f>sQcoKQe$W(>&_t$b!ErP_v9fvqF6NkKZ7#7=v66Owe+k!g&Bb{yb zxHp@wLl2DoV6j6{5iu?JfKmD{KpRgDI{jIqK;mrJ^v-Y&|tw1u_UdH9~FL;6%Kx$Q>ax>Vh*6lionMNvcxpk z=;aVRbGtx}&<|!36uo33UY6=|>-6H(AheeBFBzgUlaJ>kI@7J3<6$N zFRM46u}M)WagFD>!)i>XVzD_Nqc8e)0cG>tU*IBJ5r1aDYa>akb@NIQO4~WzDKRgJ zBOLJQcE@9kLg6Teg7?AC9;?q^GfRhsc&E2ouBZo2kRoQ#Ff=xh!QlH@=g!F&ALH&b=LzL>p%2IVyw+z0!M#;1mpRQlH*NoMI{mcJ zz(Bdx7alkC!?qSr0v;n)r5Gp_cuCApd&i#FSgGQup1(fgrohLThqQch^G{|&gL(zj zB?kSl0$1O7lCrsX9Wtrh>>25c#J(r0>_W?V5cD7DYE8Budm;fbk_&??f$Yd1;p?T4&JF(F z#URK)_7ntwLLvbVdJoJPK>Eyl_M>{Vui;x6%|o1nw*`yqPubKgVowdyay*3K;6A=X zbW%~z`V8oYmPk_BAShjmj8|Mp?uIpzdMQr@rG)R^_K7&ktQ>uDS!O2g;RBZUKpioAO)EPw8_uOLuJQX<_T{OJSuA3W~m7Ba6IdRmn zGFKF@5Z&@&qjp9b*0@8i9#T4go&)lJc;ewML98dlJL)QfN=^Jjs6|duI%Uy|?zK2& z*b1!*56uvu86V&7zCZq6QsPr$d5W*ez)*5t#Ca;W>msN^dNS>tBUqX342tz5H9vkF z1B(=hC<$TV7aM+sIwV+`^k%^^vyH5x;$atQH(+KJD$nE_e)ZYR51Q*=0p zm^L+ET|MRi6U^+Z;Ltb7oz@eXs>h~%Ii44?Mn2VoaXP7WTHB1js}1ETPz2h-N#dV( zrC^*y@r2h0WV$__c29W<5O-Dwco??CR^5R+0@W#N--hrS$#`F&dmYJQWUU8hJ7GYr zF`a=$H3=zOaZzK+yMj5z3odYK7`9I@%2Z0DSQ~+T5M%64_{*KHs=MeB8|5BEnZbZ_ zi8h|r>9x*B#4k`Z^ESn5$})F|maCoYM#N9ax$sk{RxvvhF(YB%eSb}bCGYKJ7zZ%! z@_q;qZdE%RyEt%noP2wzt^IgOoSg>(Uyq!2Ub(WL1&!99ban#%oE{kmst8rJIS)_v zzFnV|kr=(dRT}rcjTCh=T-PI!EQZOtXXJc?5e$a*J*OGnpQapPD@}k(5NnmWR~<&a zzJ%u>sgW_j5k~Lh72d5rG#_0q(ghf??x*}|lnHmqj;cn(SHQjQ4ALV&j6+#X#{p}gA+z%JDDCju1-Ij%Tz=mlx$)MSBC71*1qOAUFVUgCx5-IQ6(=hQ;t zQlpHxlX{SB_?)m3DKkTUFskh+Ph_HMBF9bz%+DpMftJ20^VhGeqw%fR>cHMmnQI)z z7lBu>%2lhM^YV)Jo2v9tcR8#H5gOzrD!=lA?)AiaArTKyUdZ&RnAK^0_ZNMF_*#=8 z0~Zoq8m`wB7gt9D$yG?D*g9s@Py)xchzj$`ka@T8E@e`A_X35`1fg7;mFUqWfkwjuC2iXrghAk) zAs!9;XWwH&2#&w%C-LQiugJ*ET1kG1h_TK<^p4a6?)l`XQwNJ!1khEFi24#9z%eqF zSs?qBEdpXtrj>L*tQU3ur~*yj_wdVQDboO$REowGC1-%cja$um^t*^0-4SQ*`;`e!7JCocUDo+USS&|5}?&+$71}0U8i|Ui|N4c#h0a7 zQLRVRQp4R}@q?ya6?n0d#uHFykVY;(0h9CY%uFzuYbsaVab3hW7`kS2669-Z@00!n zSMx5e<&Sv)-4tC>TF;l!e%daH!2He6B0}jbEoafeS9S7dl1d}8b#p~$Zl(a0_OLhA zdYr_lyydZUU$ffeYIPn+M~AM z5^jw;%v_HUX+S7M|0|r4dkCI;;+}a!&}Q;85@I8t!J=a7829T2>!i3WS0DQfp`$E~ zd1x3tIckv{&D`XU1*8H}UJMH+^(+{QgO5&4A>Y72=yb9zj2&9qgzsjN}PDt#Y30y8#C>!Ob?%{f|i8tCFHah7NyYasuaAyVYKSl7{{ zj%;w5!e*OM!OWiGXd|1C=lf-I0C`4d%)#rv5kBe`(AV?3{N=ew&u@l&dVk7YS1La% zkUap1?6TByiQ?C`Tyog1we0t};93zY^i_#o$?3r1MP*!!87Ew98s6Zn;u9e8Q0g2C z>ARIQ2O@X4t-@9?ocwM>H~&=qfEnr-;f!{wS-r@13Vg*KGY9jERwyloquR7YK$J4} zk!UyHVfp$ESXCtNAHRc9RgTq4hapFZHLOV?cJ*u4~A=<9*2#Sm(AK}X0) zuUw(#*)0-shS)N}@@M*KSqfh_6g24AO~M^ftzNW#l`$YK#ZvNvnQKt5 znt;2r+H^G_D}!@4JWk%VH>ff`p64CrKUpTApi}T18VU?76&nl;>957)ZO`K7?qhH1 z_8(=nsjKg_!HpS6QuiTxZ8%SrZbgtD%bySgVTeYGSEr*_NJ!pgjhQ<^^swV4JKnNT zTz;k>^E|vn_UA?bhd9ORwFa9|NZpW)3WLxH3Rf+&JZD@g0Ipk!&U_pr2HyFJXyh6XAw?T)d1~wfOrY-WHkkj&WwXEn~F1hFIXVsSzL9nJDHKQVqND8#d_RLd?p8%@BtdLt;ozFn*pCl8-xq4W#RYe>;~m8@rc!S5 z2MQHAfZ0p)^)|QfU+u$r8Um#4nB}Y_LBtkB1ng(j*S6V6M!!vOI<1jBVYpAhk4Q$P zJnBF1>NiS|KCxF$+e=QQ)Sh!U5`6#TJG&+hzmBzKZE5`_%N@br)QL-Cu6yAeL9sQP zJ(LgA72D#=U_$~_oT?)#s?J)dN3C0jF%C@mj_P~Xj*X6=Ii%R-1{rYO<9bt8@>L_Z z!|kHX<({QzOu5|TEE!0z%XQ>r(tYD_>{Q;4Vm>qfVhvfl`Pz;ba{QL_MwN}wYVcCO5xugUk2XwSC59Gdvm&gI6l zp>p!d?1dfzDW{~#4jh2vud3bCdx)qHyn_TMKmX~qRYZhxMC!K?`==d?*Ak3GBEy&> z=W!P1NVcK#Q`2Pgs!vlBjsRvS*;2Hh_g3CTJ*Qfzt>f`U8Lzy~_b;`Vm{_bVlk- z^G8=K82)`;y%FlePw35@w^NzI2~`os=qhDon>BBytIdlmcY47T6%Ss>h((Ab?f2J@ehzH0Tz(# znCvnP3kGI`2L^`uzdT`W>S}4B;qGeVX!Q?WRDE`E+TcZh5A1u3&sjy$68?jN}j+5%Bt7yZl^``Evm-IlYeFkI#d_iym>x6_U7W6!`X4Tt;bXw z5|a>Tv8I=^Uii~-_-7vRV`10&x>%>{#f%N9uXFPGZ=AfkeuYnKrk#qNGTSi*Y}As= zMc~(1WZv&o)}nsk9A>VF}&6e{MXzo$EpZf#>wAkmb+Y8_T>bSsZDpgkg{Q6OH22C7fzdZm;*Ceiz|Hh5l)llG%mdP? za!>Jo{-&AdaD)Fe(7<}+N&+)V)&D$8!41zFH;T({ziDRQ#f zb-3f`Si~69AK^)*Cc^ae?v->|2`Ly;hQ0}^z+)x97C9RuMzcOfe(FyTEz zJm!*|Vx9+3&f_1sXGEzkA<>o?Za8OP>#O4ha;p_)KSd%iH#N1-Q>Wg2wX;J0%}6hv zea3A#*4{@Q!L*79>dgnfBAjj1w0>pa6}WN}z3rcfYeWtwFVdGoiT2?wknKB?aD#~* z%j)DdJR13~pGMocKA@V^FDLdD6cw7dm7Gb1tU(DhYw&~_H1Rzye9E;O^ovr5P)UN= z{%yTi>A^`2k;#lHtsZP#j9QVvyc(y@E`V@mYSe0zszHsm4#t4(MWuoJ$D3$kdb5ur}T9I(Zki@F7R2kZb2K~va{F7bZOI^Gp!uL^0 z$pm=76iYZLO7nYP#HazLE|z#Tm8OaIrur!OM7B0)OzCtP8`{+1m{nVDK;?Y$hk(eN zWiYsS%Up@;4ntC{RCu|JZ275g@$)rbz_AYfN zLyny;mj8U zckoqBFVWw#{}9_`il+sv zYWHu``~Q@{Z8UB=`D7ngbvflxBmQwgNcdK3di_d6uYvOGYUP)BaDI>kQG;z)GRD*E z8Xo!dt0+RIdAH*dZ4Vsh7~C9Yw6$#S^l!F?V@MX}gcD_;tjreo&HV_i7=62Z3o42$c`2G5?o0_6RyB3f)+2<_q^V{ zZ-Ih?+kh3oOYlkf63yU>=v+99EJzlC6Ds?lzs54k87p^ZP_zbdxo9lZ5+`(^;Jn*F zE#_xXA-y=wqAGeWhNFQmMOdG+|2!{BFGgJL$InUk&8{6>@3)hoD;So`K5qVZZ-!;( z5BcbFk-L!WaxcStwskS24?OQ8l*P*DrQj*-W5D+G4Hr4cPF9(DXVo;j_4D>FyZ0xP zO}8}btvHxo9*M=R5L|`~OZ6?x7lJW^!9GI>t=;DF{G!M4xz%rvq6Xq3thhW0c4VxS z4Z-m!iTp6M$!)WqSB~_4fibHa=bO%%G1#s6#OLvZD@%?S?q@SD8u)}7FkVp0KZ@Jt zn}F{Jf0j+mCl{3;$F!y#PSY%-G^&=elq!cb4F$hO`DkpyrA2LEsiML?Kbr;PetR(h zxL|UI)bPj=@aK@sTTk>D&JMJSrM{(Hqb+CYT*sfWWo@GJ)L*{K?Tjyf_6Qz69Av{3 z>yCGH2&inpozilS2`IP&e-#5kALwuQF*p%@@_*%YXfdK zl3iNARSra`byhOs&6IK!zmJYFZmATI8Mg~}Dj-(tGJ7CDA2O;xj(y5sk>aM`corqu z9*xvTIB%+hF1(uWYk1kvELJ&Lzn8GoTDNK8K6`%i;bk-LX+BIe@jc)Im%3Bnr|);4 zb+Th!W(!Vo81VnZ?$rokVq0m@vJMM5MfR@c`T+b#<}|kPz(T@nko`!H1_uL^g#=CT z(*ijleU$-G)bJgpSCs6445oOtg*kjIC0Rri_#T!GaPG^JHKMv#KdV?A*sr)dUpmWn za`+aJo%HNhKZ!4Gr*-X2@5$;x12kHTS|H9%K7d~O9bTX|{lZ~H0Tq$i-Yh#kJjPpGnz@~)3$+z<`0#Ky?vlQqCQw37_o(1R&Z7$rOR;9Czd(nu7Xz zac^SxcqCdDm(zS_g~4yQpRgA4xW6mR$A_C;DwidPL-CX5$>Hs8Q!Bp2)@v9%t$i1- zwXdtiayTLV(EDz|vDR#+cCgS3eQ2=VfSpRBh<05wPT`x|>jp-FVPqa6IlNGLL(kw} zn_iH)l|geMb;?}aoP9$*D_KtZ#E)Qa)IYlc_|78sPPJf@^J(8bd#3R-ulYWfs%Ja! z9(RlIPoXKg^a9lXmXM#wtApQw1QnPF|1CTP1H=8R_y27sNFX#tEKrsL7U93;Q~wQW z$M^@dtAIuDZ>IXcLFOP-MFzruQ#$|MKTrD4B=(9}g#RW7z`!v6{g%JcF$Kt(0t&vx&NH!-NuJo8*tTukwyitg```NCd%y0guGPJF zo!V9B)as2@ln0)X16~CW3w%`%R8M)eUig6Y!vgbEg$kWYT=-8wDYqC;BZ~D_0|vetA;2!C0_04D(A4d-zhs# zZyMb+Q3n$ss=>eNLP;Zk8mX+g*Ko>{qidj^$|y=9-@P?TNzQSQPis!AjVj{A(EZM9 zthSzmwTIaR5s~z#29BBT7oaMAm>bFRD|A>zrYzGH8SIJ`Cax8DC+f$L-@v)6e{o(p zf+SSsvvMYQ-8!#NAMrOXQCLy{@!%GagJKi=E}@ux&6BbR#S_c6x`za(c~7k zovdSSbnL*8lu2qD-p34F>lQN&G~7g*lR+I|Nm7jm%&;<39@+PvDsUE758}ZBgYX~k zvV!3e!k}LuDyEFkw6QtG0kqK(#hji}u8;?3O5MUG;=kWN%X4AmmTvT*=Tu^0{^S+;7`M0QY zOL?N4-fhIxV~v&(NqrK~(9o07%s=H6ZyLeYhdkL8ltAEj~S zMkvrC;+E_V31kxYKJ1KD0_+LT+I)|XuJ&edy2Fn))0yRuUfMlrbRn`XVN;R%e436g*jdAf@cN65H9M&n6Y;@z544EsqpK(n1IjPAG;V z@%Gml*{nRW#>+FBqbLFtojZ(pfZz$=f;r}kIx4lEee}*Xg_Nt*1l`LhiuzqlU>=T%{FA!CNr!de@$3_((6Kr=4VxL;h1uV;>+cgls4XU)7rd8!H7NSH0; zIS#r<_H@E#{K zeoInYw-s*Kq40VCesZTa`I_)}|2?KB-@|$fDE$O=p|pHTQP3ZqaBi5V(KEo8jFOa6 zA+9z92Z$WRO|;Jfs=A z;Bz<#Fbcm=Q#sST5cBZ_dzVSlWP zu#VZ?VfS%`Oc{>%8isi9c8->iP6hU3Yt?}+O2$_TdCPZntD2) z=H_~@F1z=jU9S_C7#)$aBG@%qeHFh7NbrgQNn{e_He4?*C>VzRArbK2&t<6NcE%8P zQ^u2#?@qMIqrK@&4e00??xXU669%3AWpK~o#`V~bQH-^f~Bt-fmR+A?6caJ_82oV(z zI0a?`4j(6HJR@vrn9A{kW>hpnHJeOx@#A$lJ*m%ZkV3y?q^)&rbouFOW{_gopC;nZ z#`E1GuL2yU`U?S}%ERUo6xUh)sA%G0KwC;jOe4(_71Vf}(q%O|7p zo)4pQwZi>^nJW;wUrz*@JjE3hl14R9!iez|j8rJT%57RoFkDVof{ z9s>37?=KL;VZsbe5^7Zt$`K0ZRTAZoD9hhV7jrT)PpA@B?Ub29ZXeS!nC-4%-Jqk~fASjUOTh zRc9(c+X?Q!p-K`X`TKz7P?7>-fq}uEIU*RO$fKeX>ZlMsH2ES`s3az&oI{NM`N-8j zbUMWd5`&HyR~&!gMsdJ#LF8NK!!-BjQIflooZveExK18_qmf9GBl$;x<>g8S_lGnXyW|47jKmb!M8E-61kh!UxEas!-$kHxj7Om@uU zFig-Jr@~2Q7}?xHzQCdnY)xzDh}@z4lRN`6mtS#b$+1->;P<+}r9z2Yi?hf2C@Uvv z!pvDV6S28>PyM6mar8ac`YzV&hryWMSIJ`>mKj0SUb^-vlItN`%azL-FQyNWX<_f- z&rf^JR$iS_F_`x55zg-eh+4)IsFC?)!D@8eshmZ;)^|!@4*XdgYX;hKhG)Monid(- zN#Ku@s1k}wh8K+}@>7TNAJi5_6OwaW2d46e7wNu{y}^PG(YFvPFf2VxV8Oty3I5V6 z5%B+IR&ihI$G~4w^tSVKKUC%~ZmbDl%LW`w(;&(p*4@r1x(0|={<3SjlkefB)9{|> zB>iESOuu+WUe|rBm=zgo<4KjDV&Bv&%*~M$fHI`t7OpfZ1m>l}znu2bGvB?>?IzC`QAX}t;pQ^n<$Muo zNMC3PlF~9ndYKRAEGcxj`ART&ak)<+AQB`vS zxkXBz1&dxV+D!#jgI@=Gi;9`op(fmiu)RcQa&&2SlFyOJYymwfQPeY|_#^n*Rx0mn z2#h6?`p`FLZPa3kTVP6}r#wB_n9-g=?lML@WDG%h;DY~V?Y-n`%9Opj%(a|VX}CIO zq)5f~A)I1{#pqEos=|96PAc@Lqg4Pdn$OjO1A-Cn_oY;ZH>^t(mqXGdYK0~LLN~ds z=MIH~D8ePaW6Zj~RhdbN`v5FYP~I3J%m>X^sv@8(vujq-Fc5YRP5zuAGNRX?4x$%+ z*|oP%oFW&3|BW|2r4bPkXJgf6F*hJK4e^}@%*lO~V602zv%F}ATQ_z{UF3$3lFR8k z0YHnOjF8E7mj%!CvR4EjKv34+(ls7OA7?DSx1X}t~q%Kq^ydi;seSU)d=aCLDR9Iy}eW%wAq@q%*b5+&?go&rh8O;c< z7_4xswm6fBh}g_^cV-c`%3oxp1^&%z`Lk{>GA{M_*Vhg`VkaUz_;rOo&&wwv(C)aH zZ|5c7`3RxgHPFU!3OhGvylaOv`ZZbjRW z>&Czh?eQ!}DYHVw;clk6@spfvG}8>SF&(H!w*J*GZE`ELswbklf;i6x*Kc?E;Wm&h ztBPpF4XcZT`zJ@duP*QRCp$MRO-L{AkZO)Tl=JdgC6W(wfe6dGv`B(qs16uahi!MS zp)grcy%=PMP}#=Nc9KPvA<;^p?zh?SGlYRMh>}T%T-ftO`H1th7xMA}w9V#-dx-gA z1m9U4VMK={e)o?VAT|rEMbkEVB^28)E!5wUZNYi`K+K;|iWC+>#`7C0_Ct*$w>E__ zwtPXB2^!)#P|aI_DaTPNA%}PTDg{a>=-%T%9pUiT;!KwOHlI_HRFYVmk*5bmpk@Um zha3=2<#B}CH7AKD$CNN-i6T(o8Qfinn~%so%#MapK7`Xub_o1coADKHo5!ys9xsv$ zeARuj8Z&%~PbW|kl25fNvIyC&98VO&ln8Ei$_P#*EXYzMiZ1-|nZj7}ggtMeT$VW@ z9*tHEa|+TBK?~6&)Z_$$(*jPiSlbRb|LUKJuBB(hT5{Xm9lSd&Dg`<4Y|a`_i(UF< zt`9R(>-(R9=sC^=cdB}J^<{HQej>O^tyAw4b>W`rah+_cKL*acZ2Qr$hXsy9<>Gt* zBK%;P@P@AQd2viRh3~m*@7c8&=NE&J+u@#XG2%{#5S(cDJlPMiv`x?mnuGg5n2>w& zm>>Ru1+MgD7DK!hW|bvoX!sQTf<@~$38hD&@XLn#f+1)R08)Y+S9+ObpPmbdgn{*D ztP#vL#nD2ZMiEkuUU7FFjr2#wtWK6ne z!~-YnS5+pX)TgQjvYgJO>$YH+z(k1;ue6mLoYk2LpJE2;Pt$|mDkW0qL`{F=;HQU|U`VAF!itER({D zLnS{#TN*9r2bWo?0Nu$e6S8WbUa*JN^!stTVJz};f|E`^muT2Ce1eTb>EEAP^-6bSzCv8s1Z445xyVw_P1-o;DuHyfE~<9wE+z}fI_9nN8* zZHoET^(|IPSdxWYjE6nNNDBH)?)Dkyb&>T!mK;PvP8~Gwol~}p$JT^=11-FH0P#=9 ziRX?t+7$E=+ZmsC6jnn~j?IW`F#TfSkuh!1@p!mi^;5voUU1ZEK}1H8?3-J?R3}q! z(@jZ!s=G9&uKP3HNIOI6PLa$l-NX&hU9c27{R)6J(?)xpR8`ai26`UjoMP@qIT*){ zdT&y<K7T5y)PV&q799YQ)?FSX z2s<3AEg0LFXy-V&B$cuf*I16E08245N+u3U@F=aVW2AS&|2qcVP&E*-7*wS#EAQeo#I0VJ?Xm z3b)Z6x2pK%6GB{vQ=r%(IPd}Ai2wKh5T9rmlPZiF{E0yMryBknKwg6JJ6X+-3TaM@ z`eW{1rLz9(wa_nm2^&&S*sL%qfUxMJNa94l#&VX%pNPmYR z0nf%F0(RmlY4bt za9n6C+0G8~Q%B(H+-Yky>8bGP3^3-HT3LwbpRm>;8Q8kG)S86S-PKWRA(szHRW-X?O2lRwh45IhT^VS=o0q1#1w4i58)B@baOVjhsDd~6p zyZMd?RJ{~cMFrH3YvfY^5+>21a@}!j%2)$JF&nJNFUVVFJ7fdotcP)iWbM%^0vcsOkl)KqmYUxAg4D#zzIV{) z7-*HxtCOgfq^mZ%%3VM+;Z>ez_e|jPe$V)0gf;z*Dgya$U@8HM8IEzeAm3t2a%Xi} z#kW4=Wu8*K!ki?wlK8dooItID6Wo$}UO6<`$^!a)xA_WR^L9v(Gl>-fFOpjPDFhf; z3iLlr!|@N(fPTcRlAs1(CcF@%KH~0DtX5!|&AYaqY?LSjU@D?Qo!c^97mH@oF) zVrRtahJsmHCix(Dm)9rEYru))Yq9n<<>3LxmL^G9J}>tyFH+K>Jqx0IU$EB2vfXYU zFHN-%HFBUsq2Rex>!PaZ@xIhRlZzhkRsox4c)GZ44bR(JtMJ;YVxAW z0Fin*qvGr@_f~#ILG}e;qwS)UpY+6@UFO8&Y6>!y-x}qx zUJ_=(NxKqdL2OoPhh&jQUkUz zQV0rEhifc7MOCF?Fc2m&FP2msZN5!O^4!er=cN~3l)dNOYL;IGtEl&hsZ1-hz3C`& zW!>N3mKoqrSflsvP@DfoQw&gv8*$oRqqFhuD<%lb+58T$aSlG|DYj9`OXCh!}-VB#w2}hU2C>WGDIE!O**j4%U>m!^?eg8a9J~O;X0s^I$>g6x(PcQ^;GixtL2xPf*`*rlii! z7T1Q40YUlgUYRDpAW76tAbv&|5~=+z2A-nsn0T6ce2Kh#pCZj#M#~i`;`>K?Vc}z@ zujj4v^-O!R=lA35)sx6tjsyD;^pSZw0CNNHi08Tj28-hBw*OEMG@7 zn_;-YCCF!S>zqVl>2Ira*!^(9u}xj3Jbs&0k%Lk_3 z8gz7D^JmKUL!3rtY~ZvymH0l4rkHK^CtI_efj3~<>E^iGa)tyhK+CH@AF@bi@c#5& ztFvLI5Zs)@KjBZq67Ok(?QcUQ#h3YJ$*3X`!Xp9ZS0Fm}M z+DWC~NTp}WDslia?tQJu3pPQpUl7U^(mTWbhze71HBAvT-?sWzDQDr(?Sj#%44PDF zkz~}KDh2&M{*Z?^u(0JP=o*9p@Rc9jjrxY=8s6x$=ebFkGX=apu8wa*It74tN6)?b zHO=*OvbM04DnONs`567r_jjBP{mzf8lWXqdIqK`4zRtJj(WK7LYu@{+oi!}x4xgKw zV-U`!&;6l74B^Tu0b%kLX;}4D#kTP!RfM4*^6NFYP(~k%XkVN#+L3F~#pn=(o;Jxz zk4x2g&nf|60J@F8(`vbe__~=QI;B#w8hXt*(iPW@7r0fDI(l_wpDDKvx@`_?4V#L! zx?qwi1A9DV&(_#cUXH3VIY0_W5BtN@J)_cDEfz^w+A0>hzbL@ABUcr6cY#_N5IRdK z=4e2v%v2|46gt{tLK3)gGrnO{*;`~aFW~)?-0*K8EWV+|`varqOiE|U-^5ut2|FTw z?y)%nf!*J!Dwrw6J;P~lUyK<7We<2m?!St3A|UP}zY^*8amC3x9?djI;IhXGxGqI0 z@;j}_#SQiOQX!K{#JUAY<9Np`^N9hU}w<&C?oW1_^{`GG+8b z6KePA2H+-fB~}Zxj!AIzCXd(gHdSVJ4*u@gG+}6)`Oja)*JX5h_S028Goq1yY;i*D z6YoXL(k?x)HPmlzBSdcCPB73h+Sjk~ zylFHhhvI&<@Q=g{R-WdWYhCtnWljqoYEF9dPoJ+Prn7p2l_nEbDYse!_Ye0XmMLxC zej`+uH%FNYZPokir12M{(EY8VZaz?L!0AX>mZDXq7V2h8@CmY~RL*3EwRm#IT?Z;5 zYPOY_DMT$Xll89jaJhR`cXq05ls=_n4?zd9&uSmek@e&~eZ1v^r5@ zNiu(qv;;0U%I8jfk5t_>lVE6;XK@YXO!1XmaFz6tMKSfv0`-KWq*w|CJU2-8eHr;U@V(&UOKooV~zW|)YFQxhT{&#MEHAAw^Yy3hL)Y?rIrYBBwKy~IDa%bx%w+LRxRmzkn_S;seu>UCO z-y|%fSSILcCyt6~hvYrvM#~2PnCW4o!+u2sIG&8N`lVH2#qK}5?HG&z%@^CkXTHTp zlk-L>y5=VT%U4Z^W*ol2l5muA1R93C6Uv-%-|2j za~VXefi$DwHru{YGk*HgTEni(EHTD^HS9!H4N%TOYbXtJ#Z@1DuhKbjp+ng+h*i=- zo;`g5*D;)uaJ8r$Tx7b~-DeyRLW_uy(;vgB-+tNY&{_6&xI`mS&3LUJ$=0-SwG1+? z7odS^m(QCzTi^}SE7!h6@?5W)Cg8tnCr3hFnYUFuf&K01fVo>UJ2`t;n>hU^ zsagZ9+Ws@H&Hfa=(w**4T;yWFaUlSQZjt!))6b`XVuY$ELT0*nt^eZ{nQ~b@mBVF6 zHvFZ3?l(Wb)AsN5;a3@ZX9btx_{b{YwxPI52^$p-+Q%!GYcOs?@KB^Zr3El>v99m3 z@bV)$-gNwUJWe?POf(Tq@2g6rZlhe1>hBv-ELBm3iaHV)!TTP@FF&N!SET)0{SJ|I zM3VcEdI0&azu`C2Klh6{R=n1$V6Lv3x*B=@v1Fy+zXJoKt8N=gK?!*q4ewk7F$;3e?SUf5eNm zq%pRMQ0u91H>K~@1{%KbuSKWIkJue$+W&r9%8*y?PmqMsflVx>jDuT=Y@$e~OkFno zh89_g?_dC3-?-}F?{{84RnPyGYgXFmdxk}M$SFaENVybzD$4|)s|^%z79J`)M{8gz zV&&ny1q@PVi9U)n`?N3*C$;RIWDO`<%rja>j+0csdPf#!k0^xVi%n(0mpo=JQ45Fi z!s!`Ob+7!R;scTpqkw$~GVQ3^Zzl8)hP)GGQ7U8N!)8%WMNXt&J6p!jNQQWdHC|u# zJR2};nQxud%{y6d;I0Wk@@DjE7Ggi5pM>~qX;fF1X-ImvhS#m;NUpG4N8pD^ySa3m zOjmi!f~el9^R8TE65WfgY0F4+-acE_q|oQ#74>W!oIhHiJE|H}NzIEnUSVd>|Hv64 z2c<)&ELYW$gOg)wL-nb7nyqW|t(V^f%-PWZCok(t+pGOZ9?#K$?wid|w=xu}gZIr3 zqfIZ3Nk3>D&ui*8dR9n#0Ha@#%Y`dqZW z2uj>)T4n5#@Qk_ICCn04OX_pd+FGb(g-{DMB1cT**$c3!c8h!I2u*+_&vZ|xU3OGs z=YLoTMdfl2nXE?uQ`!-P92mWN0YWc+U}&y4tryeLAx^&qj{D!%5YpLx#pH}$qD?T5 zz}}Z)Z&YUt-UJlhJWNMeOOW?uozDh{xfJ>|=Kkz&GQ0^hV>2+Fzu4C8#E4>ynOC?u zr?WqvFqck=qmm!4n#+N%S$jt z`*^f`I5hXgf}NO3L(1@t37()yAJ3=)eoLF&jm2?K`g93*hK)ZuNhKrr;|tVrq}e-G zUru}*m&e|o_%7C20RlT`dQ@V$5)AA(>=T%Ddx2zg4i3;e>dPS>gMD9C!wn(4*ru2D z$UJ@hp)?gofHD{iKa#TcWVxGDfte1KHA%B=%fEer&@lTy1c6|?^f=-J2L{H23udbn-C8iz#^em61HyBxI42FbI0V7T5}1k8Hg zReHZIMAF-sk%_>}xw9E;BvJ1?3Oj-;}i}`>mVV;F5DIj{`zJzj7 z=w4RPtC3WoVFljf(uDOE%Y>o9SI^o~2lABzP|QbBb&D5sp5nNNTqb68j_L|eO!VylkZaMV<@XmqSwz_?T`Hi^2MhK><+@I_qP=U!;jg{+ z1UxfB1PMv*xC;y?pm6N3$&dqV^jOYDo7BC5BK!oJR^MDP5tu{>*6(P+F&V>S2m_eX z<0rwtBN|+pO{SW;yhCVcuV84m1Y*!TR>Y0Tejke(QeZtvAmp5E`2+iDXlrX(h!0%S zB)^%U=OCH%?=*-b^oFmGyjDswh+*(zNLNiM!x(H~1Hk&X(7NJz9@j#?OHv^cOz0B3 zu@^=|d(fjlxedcObv;%^ULaHr3X3BnkYWjfg{l_)3Lg9`Ov3Un#{w9rdiO-D(CBN~ zkTn)-^I--$+(ca88caLb}eV?N9#j)IWfGIMMV!>KbmvID>3wp28? zVZ)@;;0scOzlP8^yu+KEQCt$PMS4Vn7=8t)622h} zU#e_(eA!AqzYXrm(wpl)}(P?r&smc#^DbQ zwCdoLNWbE$4(w6xuwBe}ZY~sXOD&H&4WX(=s0(#vtOasINPfXl#p>FJ_6K+E7dkUQ zXh)4HenW$@j(k%5UD@=^U;~iEs!;pgh(PQ?JaC3%(*pCAyp+!&qjc(%7xw$79bDT} z*w|da^J{anVi%D5D5s*vvON}K>EjGB*3>(@V9Y^}+Y)mV>gyKNJZ6}=V(0L!czBub z@jjIXlM*;<7k^Hdz0L0mw8HrWHBVtOgM3J$NoPNbKga0u4EW)nMP&Gkm2Gi|;P+~c z^|UYarYF)+_|P#I%Hf#dp^BQJkLcR&#>XaEu?wWgP0EZZeoQ)3YpnD(sqD&!cwlfc zE3KqN(6&O_T_9;yecXoqF8HGISBS1BNWrruNDMp^fRFsXZYVAJ6W>zFS7ud!bc z)z;?^_K$RPH1cRUOto!-U9&Rw0z15X;4&6;D|1;GAn7y<>6F zEq4u`zOC#Q{EaR2r6YT;cWtBIt27&265XS3z6tZlOBm+A7Ah)_(QBF1=`(uQzgcEE zIn4DXjC3sSb`$>)cj_fK-Y+SWXb_lL0M6=U^^aqz2`*9r{SsNlgd!7>>|tslo|-o| zQ2pM?@t2;iikyCp#0_FxtlOA$Nr7vwrr|SaEboOK*mKJSP`RZ6);ZQz!%vcVamxW3 zBgUxK(;rnE4K;HK15uJ`o+kGPxb)4Y({ya}SCC*6$C7IZ38876btLGKV&Rn8z@Q7q z_kMz^@44BM*9O04%35t|^`(g`$r%w60%I)8zx|!z#IWMSv39u{^ZC$go6WDqzY2mi0oR#kgtJF_mIni;1ve_S{Q$v^aATn!3LyuM&fk zPenOjh$LY$^$r41iHg@l&>PPtZJC>G2P8<1R$*2H5lJGNAcU>N3+Fwt)h!260{cq5 zNYsPx|5(%5$%c7WL=E_i?~IlFtaBC|YR;)@eV7fCyFsw2M-Y60=7OhuGs?lzp*{=0 zAT|D2_1Jo^c0#{^_ZCKJUjf*=7ZG+Tx?k?dB*oVLJ|f#SRIT&KT)4%P{clDSFIvHY2NM>b2o4)qhXy=BI8 z<`g(MdQ-JNjvskGB-VK`liy`rJrx41C2?2x%6=}9!?h$@8eFF{7xff-%wOs(j`ts% zs5gqY9-?smPdZrZCGZs%{$?hEQt|IDnAB#9vw#Clc0^bGMl$~@g)lco9UGxnIVe)r zKiVz$tU`8-FTOqQrpw;UR~P3gRLR%32vbNaL21<~AM5(*mAH6EkR3KBw%KxG9zh|O zkDQk3;8JOQQE1>aXotXoj&qeuhB$BMMl_~=n6giT3PIj!s7HF-ws2drzwUPXWs$nE z1;Y1G-gMe#-r^I34& z8fsOcZC1}wGO}cfw|074MNbwBLPKOj$XVhCY-@XrUl5mY(pktPvgWXAs>ZX$vC3?| zLgBpn3^`1tJH5%~Ys$^==&B#bO-C};;9cV@A~!P&Nk0G2(I)&p1NDCu(s#1TPh|y^KLIa`XiMs)PR9-i87CEr$g Date: Tue, 13 Feb 2024 15:09:21 +0100 Subject: [PATCH 17/26] fix: update temporary solution input loading in Conductor Updated the instruction that manages the odd behaviour while loading input files (value set to 1 are converted in True) in method __init__. To avoid that other values are converted unexpectedly, the explicit check for value == True is added in this commit. For instance this allows to avoid the conversion of variable INDUCTANCE_MODE = 2 into INDUCTANCE_MODE = 1 (because 2 is equivalent to true but is different from it). Class: Conductor modified: conductor.py --- source_code/conductor.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source_code/conductor.py b/source_code/conductor.py index 0f04c48..bcc22ce 100644 --- a/source_code/conductor.py +++ b/source_code/conductor.py @@ -184,7 +184,7 @@ def __init__(self: Self, simulation: object, sheetConductorsList: list, ICOND: i "INDUCTANCE_MODE", "ELECTRIC_SOLVER", ] - self.operations.update({key: 1 for key in keys if self.operations[key]}) + self.operations.update({key: 1 for key in keys if self.operations[key]==True}) self.operations.update( {key: 0 for key in keys if self.operations[key] == False} ) From f2d65570280b32096a2be1c4967c84540f650778 Mon Sep 17 00:00:00 2001 From: Mortara Alessandro Date: Tue, 13 Feb 2024 16:11:46 +0100 Subject: [PATCH 18/26] fix: ValueError in method get_electric_resistance ValueError: Voltage difference along superconductor and stabilizer must be the same. Changes in method solve_current_divider. Class: StrandMixedComponent modified: strand_mixed_component.py --- source_code/strand_mixed_component.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/source_code/strand_mixed_component.py b/source_code/strand_mixed_component.py index 350908b..218d64c 100644 --- a/source_code/strand_mixed_component.py +++ b/source_code/strand_mixed_component.py @@ -641,6 +641,19 @@ def solve_current_divider( maxiter=10, disp=False, ) + + # Tolerance on newton halley increased in case I_critical is very small, + # to avoid inaccuracies on the divider that could lead to potential + # differences between sc and stab + if min(critical_current)>1e-6: + # Default tollerance in optimize.newton method + tollerance = 1.48e-8 + else: + # Value found trial and error iteration + tollerance = 1e-12 + # Other possible solution for the correct tollerance + # tollerance = min(critical_current)/1000 + # Evaluate superconducting with Halley's method sc_current = optimize.newton( self.__sc_current_residual, @@ -648,6 +661,7 @@ def solve_current_divider( args=(psi, current), fprime=self.__d_sc_current_residual, fprime2=self.__d2_sc_current_residual, + tol = tollerance, maxiter=1000, ) From ef542ec95c3deb71b8cb13674a1a53690a1dbfee Mon Sep 17 00:00:00 2001 From: Mortara Alessandro Date: Wed, 3 Apr 2024 11:20:56 +0200 Subject: [PATCH 19/26] fix: equivalent resistivity for stack_component Correction of the equation estimating the equivalent electrical resistivity of the non-SC tape materials within the method stack_electrical_resistivity_not_sc. The parallel resistor equation: rho_el_eq = A_not_sc * (sum(A_i/rho_el_i))^-1 modified: stack_component.py --- source_code/stack_component.py | 33 ++++++++++++++++++++++++++++++++- 1 file changed, 32 insertions(+), 1 deletion(-) diff --git a/source_code/stack_component.py b/source_code/stack_component.py index 2fd72e8..6d8a1a1 100644 --- a/source_code/stack_component.py +++ b/source_code/stack_component.py @@ -523,6 +523,34 @@ def stack_thermal_conductivity(self, property: dict) -> np.ndarray: axis=1 ) / self.tape_thickness + # def stack_electrical_resistivity_not_sc(self, property: dict) -> np.ndarray: + # """Method that evaluates the homogenized electrical resistivity for not superconducting materials of the stack, which is the same of the tape if the tapes constituting the stack are equals to each other. Homogenization is based on the thickness of tape layers. + + # Args: + # property (dict): dictionary with material properties in nodal points or Gauss points according to the value of flag nodal in method eval_sol_comp_properties of class SolidComponent. + + # Returns: + # np.ndarray: array with homogenized electrical resistivity of not superconducting materials of the stack of tapes in Ohm*m. + # """ + # electrical_resistivity = np.zeros( + # (property["temperature"].size, self.inputs["Material_number"] - 1) + # ) + # for ii, func in enumerate(self.electrical_resistivity_function_not_sc): + # if "cu" in func.__name__: + # electrical_resistivity[:, ii] = func( + # property["temperature"], + # property["B_field"], + # self.inputs["RRR"], + # ) + # else: + # electrical_resistivity[:, ii] = func(property["temperature"]) + # if self.inputs["Material_number"] - 1 > 1: + # # Evaluate homogenized electrical resistivity of the stack: + # # rho_el_eq = A_not_sc * (sum(A_i/rho_el_i))^-1 for any i not sc + # return self.cross_section["stab_sloped"] * np.reciprocal((self.cross_section["stab_sloped"] / electrical_resistivity).sum(axis=1)) + # elif self.inputs["Material_number"] - 1 == 1: + # return electrical_resistivity.reshape(property["temperature"].size) + def stack_electrical_resistivity_not_sc(self, property: dict) -> np.ndarray: """Method that evaluates the homogenized electrical resistivity for not superconducting materials of the stack, which is the same of the tape if the tapes constituting the stack are equals to each other. Homogenization is based on the thickness of tape layers. @@ -547,7 +575,10 @@ def stack_electrical_resistivity_not_sc(self, property: dict) -> np.ndarray: if self.inputs["Material_number"] - 1 > 1: # Evaluate homogenized electrical resistivity of the stack: # rho_el_eq = A_not_sc * (sum(A_i/rho_el_i))^-1 for any i not sc - return self.cross_section["stab_sloped"] * np.reciprocal((self.cross_section["stab_sloped"] / electrical_resistivity).sum(axis=1)) + # In this case becomes rho_el_eq = t_not_sc * (sum(t_i/rho_el_i))^-1 where t is the thickness + + return self.tape_thickness_not_sc * np.reciprocal(( + self.material_thickness_not_sc / electrical_resistivity).sum(axis=1)) elif self.inputs["Material_number"] - 1 == 1: return electrical_resistivity.reshape(property["temperature"].size) From 270e1d5da432eae1856fd0e728ef14f9a9192aae Mon Sep 17 00:00:00 2001 From: Mortara Alessandro Date: Tue, 9 Apr 2024 19:17:45 +0200 Subject: [PATCH 20/26] fix: deleted previous methos for rho_not_sc The previous method stack_electrical_resistivity_not_sc that was commented during tests of the new method for the evaluation of the electrical resistivity of the materials in parallel is now deleted. modified: stack_component.py --- source_code/stack_component.py | 28 ---------------------------- 1 file changed, 28 deletions(-) diff --git a/source_code/stack_component.py b/source_code/stack_component.py index 6d8a1a1..b4f7643 100644 --- a/source_code/stack_component.py +++ b/source_code/stack_component.py @@ -523,34 +523,6 @@ def stack_thermal_conductivity(self, property: dict) -> np.ndarray: axis=1 ) / self.tape_thickness - # def stack_electrical_resistivity_not_sc(self, property: dict) -> np.ndarray: - # """Method that evaluates the homogenized electrical resistivity for not superconducting materials of the stack, which is the same of the tape if the tapes constituting the stack are equals to each other. Homogenization is based on the thickness of tape layers. - - # Args: - # property (dict): dictionary with material properties in nodal points or Gauss points according to the value of flag nodal in method eval_sol_comp_properties of class SolidComponent. - - # Returns: - # np.ndarray: array with homogenized electrical resistivity of not superconducting materials of the stack of tapes in Ohm*m. - # """ - # electrical_resistivity = np.zeros( - # (property["temperature"].size, self.inputs["Material_number"] - 1) - # ) - # for ii, func in enumerate(self.electrical_resistivity_function_not_sc): - # if "cu" in func.__name__: - # electrical_resistivity[:, ii] = func( - # property["temperature"], - # property["B_field"], - # self.inputs["RRR"], - # ) - # else: - # electrical_resistivity[:, ii] = func(property["temperature"]) - # if self.inputs["Material_number"] - 1 > 1: - # # Evaluate homogenized electrical resistivity of the stack: - # # rho_el_eq = A_not_sc * (sum(A_i/rho_el_i))^-1 for any i not sc - # return self.cross_section["stab_sloped"] * np.reciprocal((self.cross_section["stab_sloped"] / electrical_resistivity).sum(axis=1)) - # elif self.inputs["Material_number"] - 1 == 1: - # return electrical_resistivity.reshape(property["temperature"].size) - def stack_electrical_resistivity_not_sc(self, property: dict) -> np.ndarray: """Method that evaluates the homogenized electrical resistivity for not superconducting materials of the stack, which is the same of the tape if the tapes constituting the stack are equals to each other. Homogenization is based on the thickness of tape layers. From 7feab75831aff2b1e005713ceb1fe8b42384e3ed Mon Sep 17 00:00:00 2001 From: Mortara Alessandro Date: Tue, 9 Apr 2024 21:18:22 +0200 Subject: [PATCH 21/26] fix: added StackComponent in rho_not_sc update if isinstance(strand,StrandMixedComponent): strand.dict_node_pt["electrical_resistivity_stabilizer"] = strand.strand_electrical_resistivity_not_sc( strand.dict_node_pt ) has been modified in if isinstance(strand,StrandMixedComponent) or isinstance(strand,StackComponent): strand.dict_node_pt["electrical_resistivity_stabilizer"] = strand.strand_electrical_resistivity_not_sc( strand.dict_node_pt ) to correctly update the electrical resistivity of the stabilizer in case SC is modelled as a Stack and not as a StrandMixed modified: conductor.py --- source_code/conductor.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/source_code/conductor.py b/source_code/conductor.py index a036ce9..b2866b5 100644 --- a/source_code/conductor.py +++ b/source_code/conductor.py @@ -4321,6 +4321,12 @@ def electric_solution_reorganization(self): ii :: self.inventory["StrandComponent"].number ] + obj.dict_node_pt["current_along"] = np.interp( + self.grid_features["zcoord"], + self.grid_features["zcoord_gauss"], + obj.dict_Gauss_pt["current_along"] + ) + obj.dict_Gauss_pt["delta_voltage_along"] = delta_voltage_along[ ii :: self.inventory["StrandComponent"].number ] @@ -4881,7 +4887,7 @@ def operating_conditions_em(self): else: # Update only electrical resistivity (stabilizer) at each # electric time step. - if isinstance(strand,StrandMixedComponent): + if isinstance(strand,StrandMixedComponent) or isinstance(strand,StackComponent): strand.dict_node_pt["electrical_resistivity_stabilizer"] = strand.strand_electrical_resistivity_not_sc( strand.dict_node_pt ) From 1308d1232adbc15df4c1e1fd8625c65240385665 Mon Sep 17 00:00:00 2001 From: Mortara Alessandro Date: Wed, 10 Apr 2024 10:06:29 +0200 Subject: [PATCH 22/26] fix: assignment of current_along at t=0 Assignment of current_along at t=0, considered initially set equal to op_current for each I0_OP_MODE. The variable current_along represents the current that actually flows in the strand_mix or stack component. modified: solid_component.py --- source_code/solid_component.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/source_code/solid_component.py b/source_code/solid_component.py index 07bbdb2..f33a086 100644 --- a/source_code/solid_component.py +++ b/source_code/solid_component.py @@ -267,6 +267,10 @@ def get_current(self, conductor: object): self.dict_Gauss_pt["op_current_sc"] = self.dict_Gauss_pt[ "op_current" ] + if conductor.cond_time[-1] == 0: + self.dict_node_pt["current_along"] = self.dict_node_pt["op_current"] + self.dict_Gauss_pt["current_along"] = self.dict_Gauss_pt["op_current"] + if self.flagSpecfield_current == 2: # Add also a logger @@ -297,6 +301,9 @@ def get_current(self, conductor: object): self.dict_node_pt["op_current_sc"][:-1] + self.dict_node_pt["op_current_sc"][1:] ) / 2.0 + if conductor.cond_time[-1] == 0: + self.dict_node_pt["current_along"] = self.dict_node_pt["op_current"] + self.dict_Gauss_pt["current_along"] = self.dict_Gauss_pt["op_current"] elif conductor.inputs["I0_OP_MODE"] == IOP_NOT_DEFINED: # User does not specify a current: set current carrient From 5b192ee299f1b677761c9a468ed5c3f4c85722a9 Mon Sep 17 00:00:00 2001 From: Mortara Alessandro Date: Wed, 10 Apr 2024 10:12:04 +0200 Subject: [PATCH 23/26] fix: T_cs evaluated from current_along In the previous bersion the jop was evaluated as: jop = ( np.abs(dict_dummy["op_current"]) / (self.cross_section["sc"]) ) where op_current is the current assigned by the input files. This has been changed to: jop = ( np.abs(dict_dummy["current_along"]) / (self.cross_section["sc"]) ) where current along is the fraction of current flowing inside the stack or strand component. For the calculation of T_cs, the current density was evaluated in this way. modified: strand_component.py --- source_code/strand_component.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source_code/strand_component.py b/source_code/strand_component.py index f1314a2..db4e685 100644 --- a/source_code/strand_component.py +++ b/source_code/strand_component.py @@ -270,7 +270,7 @@ def eval_tcs(self, dict_dummy): # evaluation of the current densities the total superconducting # perpendicular cross section is always used. jop = ( - np.abs(dict_dummy["op_current"]) + np.abs(dict_dummy["current_along"]) / (self.cross_section["sc"]) ) From fd78989260262d666e070f02d0973b67add27749 Mon Sep 17 00:00:00 2001 From: Mortara Alessandro Date: Wed, 10 Apr 2024 10:31:02 +0200 Subject: [PATCH 24/26] fix: Tolerance on newton halley increased Tolerance on newton halley increased in case I_critical is very small, to avoid inaccuracies on the divider that could lead to potential differences between the voltage of sc and stab. method: solve_current_divider modified: stack_component.py --- source_code/stack_component.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/source_code/stack_component.py b/source_code/stack_component.py index b4f7643..71ca45f 100644 --- a/source_code/stack_component.py +++ b/source_code/stack_component.py @@ -657,6 +657,18 @@ def solve_current_divider( disp=False, ) + # Tolerance on newton halley increased in case I_critical is very small, + # to avoid inaccuracies on the divider that could lead to potential + # differences between sc and stab + if min(critical_current)>1e-6: + # Default tollerance in optimize.newton method + tollerance = 1.48e-8 + else: + # Value found trial and error iteration + # tollerance = 1e-12 + # Other possible solution for the correct tollerance + tollerance = min(critical_current)/1000 + # Evaluate superconducting with Halley's method sc_current = optimize.newton( self.__sc_current_residual, @@ -664,6 +676,7 @@ def solve_current_divider( args=(psi, current), fprime=self.__d_sc_current_residual, fprime2=self.__d2_sc_current_residual, + tol = tollerance, maxiter=1000, ) From 4cdc096c43ee207b3379469c8c7ebf24fb5f55fa Mon Sep 17 00:00:00 2001 From: Mortara Alessandro Date: Tue, 16 Apr 2024 11:42:18 +0200 Subject: [PATCH 25/26] fix: update of electrical resistivity for Stack Inside method operating_conditions_em added: elif isinstance(strand, StackComponent): strand.dict_Gauss_pt["electrical_resistivity_stabilizer"] = strand.stack_electrical_resistivity_not_sc( strand.dict_Gauss_pt ) to update only electrical resistivity (stabilizer) at each electric time step also in presence of the StackComponent. modified: conductor.py --- source_code/conductor.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/source_code/conductor.py b/source_code/conductor.py index b2866b5..0f46f07 100644 --- a/source_code/conductor.py +++ b/source_code/conductor.py @@ -4887,10 +4887,14 @@ def operating_conditions_em(self): else: # Update only electrical resistivity (stabilizer) at each # electric time step. - if isinstance(strand,StrandMixedComponent) or isinstance(strand,StackComponent): + if isinstance(strand,StrandMixedComponent): strand.dict_node_pt["electrical_resistivity_stabilizer"] = strand.strand_electrical_resistivity_not_sc( strand.dict_node_pt ) + elif isinstance(strand, StackComponent): + strand.dict_Gauss_pt["electrical_resistivity_stabilizer"] = strand.stack_electrical_resistivity_not_sc( + strand.dict_Gauss_pt + ) elif isinstance(strand, StrandStabilizerComponent): strand.dict_node_pt["electrical_resistivity_stabilizer"] = strand.strand_electrical_resistivity( strand.dict_node_pt @@ -4979,7 +4983,11 @@ def __eval_gauss_point_em(self): # electric time step. if isinstance(strand,StrandMixedComponent): strand.dict_Gauss_pt["electrical_resistivity_stabilizer"] = strand.strand_electrical_resistivity_not_sc( - strand.dict_Gauss_pt + strand.dict_Gauss_pt + ) + elif isinstance(strand, StackComponent): + strand.dict_Gauss_pt["electrical_resistivity_stabilizer"] = strand.stack_electrical_resistivity_not_sc( + strand.dict_Gauss_pt ) elif isinstance(strand, StrandStabilizerComponent): strand.dict_Gauss_pt["electrical_resistivity_stabilizer"] = strand.strand_electrical_resistivity( From 0e3533b6edf0b2489978125b9823626e08cb4b92 Mon Sep 17 00:00:00 2001 From: Mortara Alessandro Date: Fri, 19 Apr 2024 19:33:36 +0200 Subject: [PATCH 26/26] doc: update comment in solve_current_divider Inside method solve_current_divider. Class: StackComponent modified: stack_component.py --- source_code/stack_component.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/source_code/stack_component.py b/source_code/stack_component.py index 71ca45f..1913b16 100644 --- a/source_code/stack_component.py +++ b/source_code/stack_component.py @@ -658,8 +658,9 @@ def solve_current_divider( ) # Tolerance on newton halley increased in case I_critical is very small, - # to avoid inaccuracies on the divider that could lead to potential - # differences between sc and stab + # to avoid inaccuracies on the divider that could lead to voltage + # differences between sc and stab that are not expected in the parallel + # of electric resistances. if min(critical_current)>1e-6: # Default tollerance in optimize.newton method tollerance = 1.48e-8 @@ -667,7 +668,7 @@ def solve_current_divider( # Value found trial and error iteration # tollerance = 1e-12 # Other possible solution for the correct tollerance - tollerance = min(critical_current)/1000 + tollerance = min(critical_current)/1e3 # Evaluate superconducting with Halley's method sc_current = optimize.newton(