diff --git a/examples/examples.ipynb b/examples/examples.ipynb index 8f1a21aaad..7eb8af9039 100644 --- a/examples/examples.ipynb +++ b/examples/examples.ipynb @@ -33,7 +33,7 @@ "PROJ_DIR = Path.cwd().parent\n", "\n", "\n", - "def copy_to_temp_dir(input_rel):\n", + "def copy_to_temp_dir(input_rel, config=None):\n", " \"\"\"Copy an input file to a new temp dir and return its new path.\n", "\n", " The new TemporaryDirectory object is returned to avoid destruction of the\n", @@ -61,6 +61,19 @@ " copy(input_abs_path, temp_dir_path)\n", " temp_input_path = temp_dir_path / input_abs_path.name\n", "\n", + " if config is not None:\n", + " config_rel_path = Path(config)\n", + " config_abs_path = PROJ_DIR / config_rel_path\n", + " try:\n", + " assert config_abs_path.exists()\n", + " except AssertionError as err:\n", + " raise FileNotFoundError(\"Config file doesn't exist.\") from err\n", + "\n", + " copy(config_abs_path, temp_dir_path)\n", + " temp_config_path = temp_dir_path / config_abs_path.name\n", + "\n", + " return temp_dir, temp_input_path, temp_config_path\n", + "\n", " return temp_dir, temp_input_path" ] }, @@ -86,10 +99,19 @@ "\n", "# Define input file name relative to project dir, then copy to temp dir\n", "script_dir = Path(\"__file__\").parent.resolve()\n", + "# standard run of conventional tokamak from data files\n", "input_rel = script_dir / \"data/large_tokamak_IN.DAT\"\n", - "\n", "temp_dir, temp_input_path = copy_to_temp_dir(input_rel)\n", "\n", + "# # alternative run of ST from regression files\n", + "# input_rel = script_dir / \"../tests/regression/input_files/st_regression.IN.DAT\"\n", + "# temp_dir, temp_input_path = copy_to_temp_dir(input_rel)\n", + "\n", + "# # alternative run of stellarator from regression files with associated config\n", + "# input_rel = script_dir / \"../tests/regression/input_files/stellarator_helias_eval.IN.DAT\"\n", + "# config = script_dir / \"../tests/regression/input_files/stellarator_helias_eval.stella_conf.json\"\n", + "# temp_dir, temp_input_path, temp_config = copy_to_temp_dir(input_rel, config)\n", + "\n", "# Run process on an input file in a temporary directory\n", "single_run = SingleRun(str(temp_input_path))\n", "single_run.run()" @@ -253,7 +275,7 @@ ], "metadata": { "kernelspec": { - "display_name": "env", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -270,5 +292,5 @@ } }, "nbformat": 4, - "nbformat_minor": 2 + "nbformat_minor": 4 } diff --git a/process/build.py b/process/build.py index 0de30cacc3..c939a4288f 100644 --- a/process/build.py +++ b/process/build.py @@ -802,7 +802,6 @@ def calculate_vertical_build(self, output: bool) -> None: self.outfile, "\n*Cryostat roof allowance includes uppermost PF coil and outer thermal shield.\n*Cryostat floor allowance includes lowermost PF coil, outer thermal shield and gravity support.", ) - # Output the cdivertor geometry divht = self.divgeom(output) # Issue #481 Remove build_variables.vgaptf @@ -864,6 +863,14 @@ def divgeom(self, output: bool): TART option: Peng SOFT paper """ if physics_variables.itart == 1: + if output: + po.ovarrf( + self.outfile, + "TF coil vertical offset (m)", + "(dz_tf_plasma_centre_offset)", + build_variables.dz_tf_plasma_centre_offset, + "OP ", + ) return 1.75e0 * physics_variables.rminor # Conventional tokamak divertor model # options for seperate upper and lower physics_variables.triangularity @@ -985,7 +992,6 @@ def divgeom(self, output: bool): ) divht = max(zplti, zplto) - min(zplbo, zplbi) - if output: if physics_variables.n_divertors == 1: po.oheadr(self.outfile, "Divertor build and plasma position") diff --git a/process/geometry/pfcoil_geometry.py b/process/geometry/pfcoil_geometry.py index 2b993725f0..83478b07c2 100644 --- a/process/geometry/pfcoil_geometry.py +++ b/process/geometry/pfcoil_geometry.py @@ -12,9 +12,6 @@ def pfcoil_geometry( coils_z: list[float], coils_dr: list[float], coils_dz: list[float], - dr_bore: float, - dr_cs: float, - ohdz: float, ) -> tuple[np.ndarray, np.ndarray, RectangleGeometry]: """Calculates radial and vertical distances for the geometry of the pf coils and central coil @@ -26,12 +23,6 @@ def pfcoil_geometry( :type coils_dr: List[float] :param coils_dz: list of pf coil vertical thicknesses :type coils_dz: List[float] - :param dr_bore: central solenoid inboard radius - :type dr_bore: float - :param dr_cs: central solenoid thickness - :type dr_cs: float - :param ohdz: central solenoid vertical thickness - :type ohdz: float :return: tuple containing radial and vertical coordinates for pf coils, and dataclass returning coordinates representing a rectangular geometry used to plot the central coil :rtype: Tuple[np.ndarray, np.ndarray, RectangleGeometry] """ @@ -44,9 +35,25 @@ def pfcoil_geometry( r_2 = float(coils_r[i]) + 0.5 * float(coils_dr[i]) r_points.append([r_1, r_1, r_2, r_2, r_1]) z_points.append([z_1, z_2, z_2, z_1, z_1]) + return r_points, z_points + + +def cs_geometry( + dr_bore: float, + dr_cs: float, + ohdz: float, +) -> RectangleGeometry: + """Calculates radial and vertical distances for the geometry of the central coil - central_coil = RectangleGeometry( + :param dr_bore: central solenoid inboard radius + :type dr_bore: float + :param dr_cs: central solenoid thickness + :type dr_cs: float + :param ohdz: central solenoid vertical thickness + :type ohdz: float + :return: Dataclass returning coordinates representing a rectangular geometry used to plot the central coil + :rtype: RectangleGeometry + """ + return RectangleGeometry( anchor_x=dr_bore, anchor_z=(-ohdz / 2), width=dr_cs, height=ohdz ) - - return r_points, z_points, central_coil diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 22612945c9..d2952af690 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -46,7 +46,10 @@ first_wall_geometry_double_null, first_wall_geometry_single_null, ) -from process.geometry.pfcoil_geometry import pfcoil_geometry +from process.geometry.pfcoil_geometry import ( + cs_geometry, + pfcoil_geometry, +) from process.geometry.plasma_geometry import plasma_geometry from process.geometry.shield_geometry import ( shield_geometry_double_null, @@ -2709,7 +2712,7 @@ def plot_main_plasma_information( # Add divertor information textstr_div = ( f"\n$P_{{\\text{{sep}}}}$: {mfile_data.data['p_plasma_separatrix_mw'].get_scan(scan):.2f} MW \n" - f"$\\frac{{P_{{\\text{{sep}}}}}}{{R}}$: {mfile_data.data['p_plasma_separatrix_mw/rmajor'].get_scan(scan):.2f} MW/m \n" + f"$\\frac{{P_{{\\text{{sep}}}}}}{{R}}$: {mfile_data.data['pdiv_over_r'].get_scan(scan):.2f} MW/m \n" f"$\\frac{{P_{{\\text{{sep}}}}}}{{B_T q_a R}}$: {mfile_data.data['pdivtbt_over_qar'].get_scan(scan):.2f} MW T/m " ) @@ -6992,23 +6995,39 @@ def _pack_strands_rectangular_with_obstacles( "linewidth": 2, }, ) - - textstr_superconductor = ( - f"$\\mathbf{{Superconductor:}}$\n \n" - f"Superconductor used: {sctf.SUPERCONDUCTING_TF_TYPES[mfile_data.data['i_tf_sc_mat'].get_scan(scan)]}\n" - f"Critical field at zero \ntemperature and strain: {mfile_data.data['b_tf_superconductor_critical_zero_temp_strain'].get_scan(scan):.4f} T\n" - f"Critical temperature at \nzero field and strain: {mfile_data.data['temp_tf_superconductor_critical_zero_field_strain'].get_scan(scan):.4f} K\n" - f"Temperature at conductor: {mfile_data.data['tftmp'].get_scan(scan):.4f} K\n" - f"$I_{{\\text{{TF,turn critical}}}}$: {mfile_data.data['c_turn_cables_critical'].get_scan(scan):,.2f} A\n" - f"$I_{{\\text{{TF,turn}}}}$: {mfile_data.data['c_tf_turn'].get_scan(scan):,.2f} A\n" - f"Critcal current ratio: {mfile_data.data['f_c_tf_turn_operating_critical'].get_scan(scan):,.4f}\n" - f"Superconductor temperature \nmargin: {mfile_data.data['temp_tf_superconductor_margin'].get_scan(scan):,.4f} K\n" - f"\n$\\mathbf{{Quench:}}$\n \n" - f"Quench dump time: {mfile_data.data['t_tf_superconductor_quench'].get_scan(scan):.4e} s\n" - f"Quench detection time: {mfile_data.data['t_tf_quench_detection'].get_scan(scan):.4e} s\n" - f"User input max temperature \nduring quench: {mfile_data.data['temp_tf_conductor_quench_max'].get_scan(scan):.2f} K\n" - f"Required maxium WP current \ndensity for heat protection:\n{mfile_data.data['j_tf_wp_quench_heat_max'].get_scan(scan):.2e} A/m$^2$\n" - ) + if mfile_data.data["i_tf_sc_mat"].get_scan(scan) > 5: + textstr_superconductor = ( + f"$\\mathbf{{Superconductor:}}$\n \n" + f"Superconductor used: {sctf.SUPERCONDUCTING_TF_TYPES[mfile_data.data['i_tf_sc_mat'].get_scan(scan)]}\n" + f"Critical field at zero \ntemperature and strain: {mfile_data.data['b_tf_superconductor_critical_zero_temp_strain'].get_scan(scan):.4f} T\n" + f"Critical temperature at \nzero field and strain: {mfile_data.data['temp_tf_superconductor_critical_zero_field_strain'].get_scan(scan):.4f} K\n" + f"Temperature at conductor: {mfile_data.data['tftmp'].get_scan(scan):.4f} K\n" + f"$I_{{\\text{{TF,turn critical}}}}$: {mfile_data.data['c_turn_cables_critical'].get_scan(scan):,.2f} A\n" + f"$I_{{\\text{{TF,turn}}}}$: {mfile_data.data['c_tf_turn'].get_scan(scan):,.2f} A\n" + f"Critcal current ratio: {mfile_data.data['f_c_tf_turn_operating_critical'].get_scan(scan):,.4f}\n" + f"Superconductor temperature \nmargin: {mfile_data.data['temp_tf_superconductor_margin'].get_scan(scan):,.4f} K\n" + f"\n$\\mathbf{{Quench:}}$\n \n" + f"Quench dump time: {mfile_data.data['t_tf_superconductor_quench'].get_scan(scan):.4e} s\n" + f"Quench detection time: {mfile_data.data['t_tf_quench_detection'].get_scan(scan):.4e} s\n" + f"Required maxium WP current \ndensity for heat protection:\n{mfile_data.data['j_tf_wp_quench_heat_max'].get_scan(scan):.2e} A/m$^2$\n" + ) + else: + textstr_superconductor = ( + f"$\\mathbf{{Superconductor:}}$\n \n" + f"Superconductor used: {sctf.SUPERCONDUCTING_TF_TYPES[mfile_data.data['i_tf_sc_mat'].get_scan(scan)]}\n" + f"Critical field at zero \ntemperature and strain: {mfile_data.data['b_tf_superconductor_critical_zero_temp_strain'].get_scan(scan):.4f} T\n" + f"Critical temperature at \nzero field and strain: {mfile_data.data['temp_tf_superconductor_critical_zero_field_strain'].get_scan(scan):.4f} K\n" + f"Temperature at conductor: {mfile_data.data['tftmp'].get_scan(scan):.4f} K\n" + f"$I_{{\\text{{TF,turn critical}}}}$: {mfile_data.data['c_turn_cables_critical'].get_scan(scan):,.2f} A\n" + f"$I_{{\\text{{TF,turn}}}}$: {mfile_data.data['c_tf_turn'].get_scan(scan):,.2f} A\n" + f"Critcal current ratio: {mfile_data.data['f_c_tf_turn_operating_critical'].get_scan(scan):,.4f}\n" + f"Superconductor temperature \nmargin: {mfile_data.data['temp_tf_superconductor_margin'].get_scan(scan):,.4f} K\n" + f"\n$\\mathbf{{Quench:}}$\n \n" + f"Quench dump time: {mfile_data.data['t_tf_superconductor_quench'].get_scan(scan):.4e} s\n" + f"Quench detection time: {mfile_data.data['t_tf_quench_detection'].get_scan(scan):.4e} s\n" + f"User input max temperature \nduring quench: {mfile_data.data['temp_tf_conductor_quench_max'].get_scan(scan):.2f} K\n" + f"Required maxium WP current \ndensity for heat protection:\n{mfile_data.data['j_tf_wp_quench_heat_max'].get_scan(scan):.2e} A/m$^2$\n" + ) axis.text( 0.75, 0.9, @@ -7130,9 +7149,13 @@ def plot_pf_coils(axis, mfile_data, scan, colour_scheme): coils_dz = [] coil_text = [] + # Check for Central Solenoid + iohcl = mfile_data.data["iohcl"].get_scan(scan) if "iohcl" in mfile_data.data else 1 + dr_bore = mfile_data.data["dr_bore"].get_scan(scan) - dr_cs = mfile_data.data["dr_cs"].get_scan(scan) - dz_cs_full = mfile_data.data["dz_cs_full"].get_scan(scan) + if iohcl == 1: + dr_cs = mfile_data.data["dr_cs"].get_scan(scan) + dz_cs_full = mfile_data.data["dz_cs_full"].get_scan(scan) # Number of coils, both PF and CS number_of_coils = 0 @@ -7140,9 +7163,6 @@ def plot_pf_coils(axis, mfile_data, scan, colour_scheme): if "r_pf_coil_middle[" in item: number_of_coils += 1 - # Check for Central Solenoid - iohcl = mfile_data.data["iohcl"].get_scan(scan) if "iohcl" in mfile_data.data else 1 - # If Central Solenoid present, ignore last entry in for loop # The last entry will be the OH coil in this case noc = number_of_coils - 1 if iohcl == 1 else number_of_coils @@ -7154,28 +7174,32 @@ def plot_pf_coils(axis, mfile_data, scan, colour_scheme): coils_dz.append(mfile_data.data[f"pfdz({coil:01})"].get_scan(scan)) coil_text.append(str(coil + 1)) - r_points, z_points, central_coil = pfcoil_geometry( + r_points, z_points = pfcoil_geometry( coils_r=coils_r, coils_z=coils_z, coils_dr=coils_dr, coils_dz=coils_dz, - dr_bore=dr_bore, - dr_cs=dr_cs, - ohdz=dz_cs_full, ) - # Plot CS compression structure - r_precomp_outer, r_precomp_inner = cumulative_radial_build2( - "dr_cs_precomp", mfile_data, scan - ) - axis.add_patch( - patches.Rectangle( - xy=(r_precomp_inner, central_coil.anchor_z), - width=r_precomp_outer - r_precomp_inner, - height=central_coil.height, - facecolor=CSCOMPRESSION_COLOUR[colour_scheme - 1], + if iohcl == 1: + central_coil = cs_geometry( + dr_bore=dr_bore, + dr_cs=dr_cs, + ohdz=dz_cs_full, + ) + + # Plot CS compression structure + r_precomp_outer, r_precomp_inner = cumulative_radial_build2( + "dr_cs_precomp", mfile_data, scan + ) + axis.add_patch( + patches.Rectangle( + xy=(r_precomp_inner, central_coil.anchor_z), + width=r_precomp_outer - r_precomp_inner, + height=central_coil.height, + facecolor=CSCOMPRESSION_COLOUR[colour_scheme - 1], + ) ) - ) # Get axis height for fontsize scaling bbox = axis.get_window_extent().transformed(axis.figure.dpi_scale_trans.inverted()) @@ -7193,14 +7217,15 @@ def plot_pf_coils(axis, mfile_data, scan, colour_scheme): va="center", fontsize=fontsize, ) - axis.add_patch( - patches.Rectangle( - xy=(central_coil.anchor_x, central_coil.anchor_z), - width=central_coil.width, - height=central_coil.height, - facecolor=SOLENOID_COLOUR[colour_scheme - 1], + if iohcl == 1: + axis.add_patch( + patches.Rectangle( + xy=(central_coil.anchor_x, central_coil.anchor_z), + width=central_coil.width, + height=central_coil.height, + facecolor=SOLENOID_COLOUR[colour_scheme - 1], + ) ) - ) def plot_info(axis, data, mfile_data, scan): @@ -7595,24 +7620,49 @@ def plot_magnetics_info(axis, mfile_data, scan): ].get_scan(scan) if i_tf_sup == 1: - data = [ - (pf_info[0][0], pf_info[0][1], "MA"), - (pf_info[1][0], pf_info[1][1], "MA"), - (pf_info_3_a, pf_info_3_b, "MA"), - (vssoft, "Startup flux swing", "Wb"), - ("vs_cs_pf_total_pulse", "Available flux swing", "Wb"), - (t_plant_pulse_burn, "Burn time", "hrs"), - ("", "", ""), - (f"#TF coil type is {tftype}", "", ""), - ("b_tf_inboard_peak_with_ripple", "Peak field at conductor (w. rip.)", "T"), - ("f_c_tf_turn_operating_critical", r"I/I$_{\mathrm{crit}}$", ""), - ("temp_tf_superconductor_margin", "TF Temperature margin", "K"), - ("temp_cs_superconductor_margin", "CS Temperature margin", "K"), - (sig_cond, "TF Cond max TRESCA stress", "MPa"), - (sig_case, "TF Case max TRESCA stress", "MPa"), - ("m_tf_coils_total/n_tf_coils", "Mass per TF coil", "kg"), - ] - + if mfile_data.data["iohcl"].get_scan(scan) == 1: + data = [ + (pf_info[0][0], pf_info[0][1], "MA"), + (pf_info[1][0], pf_info[1][1], "MA"), + (pf_info_3_a, pf_info_3_b, "MA"), + (vssoft, "Startup flux swing", "Wb"), + ("vs_cs_pf_total_pulse", "Available flux swing", "Wb"), + (t_plant_pulse_burn, "Burn time", "hrs"), + ("", "", ""), + (f"#TF coil type is {tftype}", "", ""), + ( + "b_tf_inboard_peak_with_ripple", + "Peak field at conductor (w. rip.)", + "T", + ), + ("f_c_tf_turn_operating_critical", r"I/I$_{\mathrm{crit}}$", ""), + ("temp_tf_superconductor_margin", "TF Temperature margin", "K"), + ("temp_cs_superconductor_margin", "CS Temperature margin", "K"), + (sig_cond, "TF Cond max TRESCA stress", "MPa"), + (sig_case, "TF Case max TRESCA stress", "MPa"), + ("m_tf_coils_total/n_tf_coils", "Mass per TF coil", "kg"), + ] + else: + data = [ + (pf_info[0][0], pf_info[0][1], "MA"), + (pf_info[1][0], pf_info[1][1], "MA"), + (pf_info_3_a, pf_info_3_b, "MA"), + (vssoft, "Startup flux swing", "Wb"), + ("vs_cs_pf_total_pulse", "Available flux swing", "Wb"), + (t_plant_pulse_burn, "Burn time", "hrs"), + ("", "", ""), + (f"#TF coil type is {tftype}", "", ""), + ( + "b_tf_inboard_peak_with_ripple", + "Peak field at conductor (w. rip.)", + "T", + ), + ("f_c_tf_turn_operating_critical", r"I/I$_{\mathrm{crit}}$", ""), + ("temp_tf_superconductor_margin", "TF Temperature margin", "K"), + (sig_cond, "TF Cond max TRESCA stress", "MPa"), + (sig_case, "TF Case max TRESCA stress", "MPa"), + ("m_tf_coils_total/n_tf_coils", "Mass per TF coil", "kg"), + ] else: p_cp_resistive = 1.0e-6 * mfile_data.data["p_cp_resistive"].get_scan(scan) p_tf_leg_resistive = 1.0e-6 * mfile_data.data["p_tf_leg_resistive"].get_scan( @@ -12678,13 +12728,13 @@ def main_plot( plot_density_limit_comparison(fig20.add_subplot(221), m_file_data, scan) plot_confinement_time_comparison(fig20.add_subplot(224), m_file_data, scan) plot_current_profiles_over_time(fig21.add_subplot(111), m_file_data, scan) - - plot_cs_coil_structure( - fig22.add_subplot(121, aspect="equal"), fig22, m_file_data, scan - ) - plot_cs_turn_structure( - fig22.add_subplot(224, aspect="equal"), fig22, m_file_data, scan - ) + if m_file_data.data["iohcl"].get_scan(scan) == 1: + plot_cs_coil_structure( + fig22.add_subplot(121, aspect="equal"), fig22, m_file_data, scan + ) + plot_cs_turn_structure( + fig22.add_subplot(224, aspect="equal"), fig22, m_file_data, scan + ) plot_first_wall_top_down_cross_section( fig23.add_subplot(221, aspect="equal"), m_file_data, scan @@ -12967,7 +13017,7 @@ def main(args=None): if int(m_file.data["i_single_null"].get_scan(scan)) == 0: vertical_upper = [ "z_plasma_xpoint_upper", - "dz_fw_plasma_gap", + "dz_xpoint_divertor", "dz_divertor", "dz_shld_upper", "dz_vv_upper", @@ -12975,6 +13025,7 @@ def main(args=None): "dz_shld_thermal", "dr_tf_shld_gap", "dr_tf_inboard", + "dz_tf_cryostat", ] else: vertical_upper = [ diff --git a/process/physics.py b/process/physics.py index 8d16bf4559..41732fbb0d 100644 --- a/process/physics.py +++ b/process/physics.py @@ -5890,15 +5890,15 @@ def outplas(self): # Double null divertor configuration po.ovarre( self.outfile, - "Pdivt / R ratio (MW/m) (On peak divertor)", - "(p_div_separatrix_max_mw/physics_variables.rmajor)", + "p_plasma_separatrix_mw / R ratio (MW/m) (On peak divertor)", + "(pdiv_over_r)", physics_variables.p_div_separatrix_max_mw / physics_variables.rmajor, "OP ", ) po.ovarre( self.outfile, - "Pdivt Bt / qAR ratio (MWT/m) (On peak divertor)", - "(pdivmaxbt/qar)", + "p_div_separatrix_max_mw * Bt / q A R ratio (MWT/m) (On peak divertor)", + "(pdivtbt_over_qar)", ( ( physics_variables.p_div_separatrix_max_mw @@ -5916,14 +5916,14 @@ def outplas(self): # Single null divertor configuration po.ovarre( self.outfile, - "Psep / R ratio (MW/m)", - "(p_plasma_separatrix_mw/rmajor)", + "p_plasma_separatrix_mw / R ratio (MW/m)", + "(pdiv_over_r)", physics_variables.p_plasma_separatrix_mw / physics_variables.rmajor, "OP ", ) po.ovarre( self.outfile, - "Psep Bt / qAR ratio (MWT/m)", + "p_div_separatrix_max_mw * Bt / q A R ratio (MWT/m)", "(pdivtbt_over_qar)", ( (