From 5e89fecff5671ecc784a0c7ca2bf08cb5b52e58f Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Mon, 6 May 2024 18:22:52 -0500 Subject: [PATCH] Restructure NWL files --- Examples/nwl_geom_example.py | 64 ++++++ Examples/nwl_plot_example.py | 21 ++ ...example_python.py => parastell_example.py} | 0 NWL/NWL_geom.py | 35 --- NWL/NWL_plot.py | 19 -- NWL/NWL_transport.py | 11 - NWL/NWL.py => parastell/nwl_utils.py | 206 ++---------------- 7 files changed, 99 insertions(+), 257 deletions(-) create mode 100644 Examples/nwl_geom_example.py create mode 100644 Examples/nwl_plot_example.py rename Examples/{example_python.py => parastell_example.py} (100%) delete mode 100644 NWL/NWL_geom.py delete mode 100644 NWL/NWL_plot.py delete mode 100644 NWL/NWL_transport.py rename NWL/NWL.py => parastell/nwl_utils.py (52%) diff --git a/Examples/nwl_geom_example.py b/Examples/nwl_geom_example.py new file mode 100644 index 00000000..246aed21 --- /dev/null +++ b/Examples/nwl_geom_example.py @@ -0,0 +1,64 @@ +from pathlib import Path + +import numpy as np + +import parastell.parastell as ps + + +# Define directory to export all output files to +export_dir = '' +# Define plasma equilibrium VMEC file +vmec_file = 'wout_vmec.nc' + +# Instantiate ParaStell build +stellarator = ps.Stellarator(vmec_file) + +# Define build parameters for in-vessel components +toroidal_angles = [0.0, 30.0, 60.0, 90.0] +poloidal_angles = [0.0, 120.0, 240.0, 360.0] +wall_s = 1.08 + +empty_radial_build_dict = {} +# Construct in-vessel components +stellarator.construct_invessel_build( + toroidal_angles, + poloidal_angles, + wall_s, + empty_radial_build_dict, + split_chamber=False +) +# Export in-vessel component files +stellarator.export_invessel_build( + export_cad_to_dagmc=False, + export_dir=export_dir +) + +# Define source mesh parameters +mesh_size = (11, 81, 61) +toroidal_extent = 90.0 +# Construct source +stellarator.construct_source_mesh( + mesh_size, + toroidal_extent +) +# Export source file +stellarator.export_source_mesh( + filename='source_mesh', + export_dir=export_dir +) + +strengths = stellarator.source_mesh.strengths + +filepath = Path(export_dir) / 'strengths.txt' + +file = open(filepath, 'w') +for tet in strengths: + file.write(f'{tet}\n') + +# Export DAGMC neutronics H5M file +stellarator.export_dagmc( + skip_imprint=False, + legacy_faceting=True, + filename='dagmc', + export_dir=export_dir +) diff --git a/Examples/nwl_plot_example.py b/Examples/nwl_plot_example.py new file mode 100644 index 00000000..68f59f22 --- /dev/null +++ b/Examples/nwl_plot_example.py @@ -0,0 +1,21 @@ +import parastell.nwl_utils as nwl + + +# Define simulation parameters +dagmc_geom = 'first_wall.h5m' +source_mesh = 'SourceMesh.h5m' +tor_ext = 90.0 +ss_file = 'strengths.txt' +num_parts = 100000 + +nwl.nwl_transport(dagmc_geom, source_mesh, tor_ext, ss_file, num_parts) + +# Define first wall geometry and plotting parameters +source_file = 'surface_source.h5' +ss_file = 'strengths.txt' +plas_eq = 'plas_eq.nc' +tor_ext = 90.0 +pol_ext = 360.0 +wall_s = 1.08 + +nwl.nwl_plot(source_file, ss_file, plas_eq, tor_ext, pol_ext, wall_s) diff --git a/Examples/example_python.py b/Examples/parastell_example.py similarity index 100% rename from Examples/example_python.py rename to Examples/parastell_example.py diff --git a/NWL/NWL_geom.py b/NWL/NWL_geom.py deleted file mode 100644 index aaa272ec..00000000 --- a/NWL/NWL_geom.py +++ /dev/null @@ -1,35 +0,0 @@ -import NWL.NWL as NWL -import numpy as np - - -# Define first wall geometry parameters -plas_eq = 'plas_eq.nc' -wall_s = 1.2 -tor_ext = 90.0 -num_phi = 61 -num_theta = 61 - -# Define fusion neutron source parameters -source = { - 'num_s': 11, - 'num_theta': 81, - 'num_phi': 61, - 'tor_ext': tor_ext -} - -# Define export parameters -export = { - 'step_export': True, - 'h5m_export': 'Cubit', - 'h5m_filename': 'first_wall', - 'dir': '', - 'native_meshing': False, - 'facet_tol': 1, - 'len_tol': 5, - 'norm_tol': None -} - -NWL.NWL_geom( - plas_eq, wall_s, tor_ext, num_phi, num_theta, source = source, - export = export -) diff --git a/NWL/NWL_plot.py b/NWL/NWL_plot.py deleted file mode 100644 index d3fff406..00000000 --- a/NWL/NWL_plot.py +++ /dev/null @@ -1,19 +0,0 @@ -import NWL.NWL as NWL -import numpy as np - - -# Define first wall geometry and plotting parameters -source_file = 'surface_source.h5' -ss_file = 'strengths.txt' -plas_eq = 'plas_eq.nc' -tor_ext = 90.0 -pol_ext = 360.0 -num_phi = 101 -num_theta = 101 -wall_s = 1.2 -num_levels = 10 - -NWL.NWL_plot( - source_file, ss_file, plas_eq, tor_ext, pol_ext, wall_s, num_phi, - num_theta, num_levels -) diff --git a/NWL/NWL_transport.py b/NWL/NWL_transport.py deleted file mode 100644 index 6a83cc66..00000000 --- a/NWL/NWL_transport.py +++ /dev/null @@ -1,11 +0,0 @@ -import NWL.NWL as NWL - - -# Define simulation parameters -dagmc_geom = 'first_wall.h5m' -source_mesh = 'SourceMesh.h5m' -tor_ext = 90.0 -ss_file = 'strengths.txt' -num_parts = 100000 - -NWL.NWL_transport(dagmc_geom, source_mesh, tor_ext, ss_file, num_parts) diff --git a/NWL/NWL.py b/parastell/nwl_utils.py similarity index 52% rename from NWL/NWL.py rename to parastell/nwl_utils.py index 6a3d79a7..4fa8ab11 100644 --- a/NWL/NWL.py +++ b/parastell/nwl_utils.py @@ -1,177 +1,9 @@ +import h5py import numpy as np - - -# Define default export dictionary -export_def = { - 'step_export': True, - 'h5m_export': None, - 'dir': '', - 'h5m_filename': 'dagmc', - 'native_meshing': False, - 'facet_tol': None, - 'len_tol': None, - 'norm_tol': None, - 'anisotropic_ratio': 100, - 'deviation_angle': 5, - 'min_mesh_size': 5.0, - 'max_mesh_size': 20.0, - 'volume_atol': 0.00001, - 'center_atol': 0.00001, - 'bounding_box_atol': 0.00001 -} - - -def NWL_geom( - plas_eq, wall_s, tor_ext, num_phi = 61, num_theta = 61, source = None, - export = export_def, logger = None - ): - """Creates DAGMC-compatible neutronics H5M geometry of first wall. - - Arguments: - plas_eq (str): path to plasma equilibrium NetCDF file. - wall_s (float): closed flux surface label extrapolation at wall. - tor_ext (float): toroidal extent to model (deg). - num_phi (int): number of phi geometric cross-sections to make for each - build segment (defaults to 61). - num_theta (int): number of points defining the geometric cross-section - (defaults to 61). - source (dict): dictionary of source mesh parameters including - { - 'num_s': number of closed magnetic flux surfaces defining mesh - (int), - 'num_theta': number of poloidal angles defining mesh (int), - 'num_phi': number of toroidal angles defining mesh (int) - } - export (dict): dictionary of export parameters including - { - 'step_export': export component STEP files (bool, defaults to - True), - 'h5m_export': export DAGMC-compatible neutronics H5M file using - Cubit or Gmsh. Acceptable values are None or a string value - of 'Cubit' or 'Gmsh' (str, defaults to None). The string is - case-sensitive. Note that if magnets are included, 'Cubit' - must be used, - 'dir': directory to which to export output files (str, defaults - to empty string). Note that directory must end in '/', if - using Linux or MacOS, or '\' if using Windows. - 'h5m_filename': name of DAGMC-compatible neutronics H5M file - (str, defaults to 'dagmc'), - 'native_meshing': choose native or legacy faceting for DAGMC - export (bool, defaults to False), - 'facet_tol': maximum distance a facet may be from surface of - CAD representation for Cubit export (float, defaults to - None), - 'len_tol': maximum length of facet edge for Cubit export - (float, defaults to None), - 'norm_tol': maximum change in angle between normal vector of - adjacent facets (float, defaults to None), - 'anisotropic_ratio': controls edge length ratio of elements - (float, defaults to 100.0), - 'deviation_angle': controls deviation angle of facet from - surface, i.e. lower deviation angle => more elements in - areas with higher curvature (float, defaults to 5.0), - 'min_mesh_size': minimum mesh element size for Gmsh export - (float, defaults to 5.0), - 'max_mesh_size': maximum mesh element size for Gmsh export - (float, defaults to 20.0), - 'volume_atol': absolute volume tolerance to allow when matching - parts in intermediate BREP file with CadQuery parts for - Gmsh export(float, defaults to 0.00001), - 'center_atol': absolute center coordinates tolerance to allow - when matching parts in intermediate BREP file with CadQuery - parts for Gmsh export (float, defaults to 0.00001), - 'bounding_box_atol': absolute bounding box tolerance to allow - when matching parts in intermediate BREP file with CadQuery - parts for Gmsh export (float, defaults to 0.00001). - } - logger (object): logger object (defaults to None). If no logger is - supplied, a default logger will be instantiated. - - Returns: - strengths (list): list of source strengths for each tetrahedron (1/s). - Returned only if source mesh is generated. - """ - import parastell as ps - import source_mesh as sm - import parastell.log as log - import read_vmec - import cubit - from scipy.interpolate import RegularGridInterpolator - import os - import inspect - from pathlib import Path - - export_dict = export_def.copy() - export_dict.update(export) - - if export_dict['h5m_export'] == 'Cubit': - cubit_dir = os.path.dirname(inspect.getfile(cubit)) - cubit_dir = Path(cubit_dir) / Path('plugins') - cubit.init([ - 'cubit', - '-nojournal', - '-nographics', - '-information', 'off', - '-warning', 'off', - '-commandplugindir', - str(cubit_dir) - ]) - - if logger == None or not logger.hasHandlers(): - logger = log.init() - - logger.info('Building first wall geometry...') - - vmec = read_vmec.vmec_data(plas_eq) - - repeat = 0 - - # Define arrays of toroidal and poloidal angles of first wall build - tor_ext = np.deg2rad(tor_ext) - phi_list = np.linspace(0, tor_ext, num = num_phi) - theta_list = np.linspace(0, 2*np.pi, num = num_theta) - - # Define offset matrix - offset_mat = np.zeros((num_phi, num_theta)) - - # Define volume used to cut periods - cutter = None - - # Build offset interpolator - interp = RegularGridInterpolator((phi_list, theta_list), offset_mat) - - # Generate stellarator torus with exterior surface at first wall - try: - torus, cutter = ps.stellarator_torus( - vmec, wall_s, tor_ext, repeat, phi_list, theta_list, interp, cutter - ) - except ValueError as e: - logger.error(e.args[0]) - raise e - - components = { - 'first_wall': { - 'solid': torus, - 'h5m_tag': 'Vacuum' - } - } - - try: - ps.exports(export, components, None, logger) - except ValueError as e: - logger.error(e.args[0]) - raise e - - if source is not None: - strengths = sm.source_mesh(vmec, source, export['dir'], logger = logger) - - filepath = Path(export['dir']) / 'strengths.txt' - - file = open(filepath, 'w') - for tet in strengths: - file.write(f'{tet}\n') - - return strengths +from scipy.optimize import direct +import openmc +import src.pystell.read_vmec as read_vmec +import matplotlib.pyplot as plt def extract_ss(ss_file): @@ -194,7 +26,7 @@ def extract_ss(ss_file): return strengths -def NWL_transport(dagmc_geom, source_mesh, tor_ext, ss_file, num_parts): +def nwl_transport(dagmc_geom, source_mesh, tor_ext, ss_file, num_parts): """Performs neutron transport on first wall geometry via OpenMC. Arguments: @@ -204,8 +36,6 @@ def NWL_transport(dagmc_geom, source_mesh, tor_ext, ss_file, num_parts): ss_file (str): source strength input file. num_parts (int): number of source particles to simulate. """ - import openmc - tor_ext = np.deg2rad(tor_ext) strengths = extract_ss(ss_file) @@ -302,8 +132,6 @@ def find_coords(vmec, wall_s, phi, pt): Returns: theta (float): poloidal angle (rad). """ - from scipy.optimize import direct - # Solve for the poloidal angle via minimization theta = direct( min_problem, @@ -352,8 +180,6 @@ def extract_coords(source_file): coords (array of array of float): Cartesian coordinates of all particle surface crossings. """ - import h5py - # Load source file file = h5py.File(source_file, 'r') # Extract source information @@ -367,24 +193,22 @@ def extract_coords(source_file): return coords -def plot(NWL_mat, phi_pts, theta_pts, num_levels): +def plot(nwl_mat, phi_pts, theta_pts, num_levels): """Generates contour plot of NWL. Arguments: - NWL_mat (array of array of float): NWL solutions at centroids of + nwl_mat (array of array of float): NWL solutions at centroids of (phi, theta) bins (MW). phi_pts (array of float): centroids of toroidal angle bins (rad). theta_bins (array of float): centroids of poloidal angle bins (rad). num_levels (int): number of contour regions. """ - import matplotlib.pyplot as plt - phi_pts = np.rad2deg(phi_pts) theta_pts = np.rad2deg(theta_pts) - levels = np.linspace(np.min(NWL_mat), np.max(NWL_mat), num = num_levels) + levels = np.linspace(np.min(nwl_mat), np.max(nwl_mat), num = num_levels) fig, ax = plt.subplots() - CF = ax.contourf(phi_pts, theta_pts, NWL_mat, levels = levels) + CF = ax.contourf(phi_pts, theta_pts, nwl_mat, levels = levels) cbar = plt.colorbar(CF) cbar.ax.set_ylabel('NWL (MW)') plt.xlabel('Toroidal Angle (degrees)') @@ -392,7 +216,7 @@ def plot(NWL_mat, phi_pts, theta_pts, num_levels): fig.savefig('NWL.png') -def NWL_plot( +def nwl_plot( source_file, ss_file, plas_eq, tor_ext, pol_ext, wall_s, num_phi = 101, num_theta = 101, num_levels = 10 ): @@ -409,15 +233,13 @@ def NWL_plot( num_theta (int): number of poloidal angle bins (defaults to 101). num_levels (int): number of contour regions (defaults to 10). """ - import read_vmec - tor_ext = np.deg2rad(tor_ext) pol_ext = np.deg2rad(pol_ext) coords = extract_coords(source_file) # Load plasma equilibrium data - vmec = read_vmec.vmec_data(plas_eq) + vmec = read_vmec.VMECData(plas_eq) phi_coords, theta_coords = flux_coords(vmec, wall_s, coords) @@ -451,6 +273,6 @@ def NWL_plot( # Define number of source particles num_parts = len(coords) - NWL_mat = count_mat*n_energy*eV2J*SS*J2MJ/num_parts + nwl_mat = count_mat*n_energy*eV2J*SS*J2MJ/num_parts - plot(NWL_mat, phi_pts, theta_pts, num_levels) + plot(nwl_mat, phi_pts, theta_pts, num_levels)