diff --git a/pyproject.toml b/pyproject.toml index f54fa1d2..56a62598 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -46,9 +46,11 @@ dependencies = [ "pyyaml", "megpy>=2.0.3", "fibe>=0.3", + "periodictable", "vmecpp", - # "quends @ git+https://github.com/sandialabs/quends.git" "scikit-image", # Stricly not for MITIM, but good to have for pygacode + # "extended-lengyel @ git+https://github.com/cfs-energy/extended-lengyel.git", + # "quends @ git+https://github.com/sandialabs/quends.git" #"onnx2pytorch", # Stricly not for MITIM, but good to use ONNX models ] diff --git a/src/mitim_modules/maestro/MAESTROmain.py b/src/mitim_modules/maestro/MAESTROmain.py index 403f2282..bc49e2a1 100644 --- a/src/mitim_modules/maestro/MAESTROmain.py +++ b/src/mitim_modules/maestro/MAESTROmain.py @@ -9,9 +9,10 @@ from mitim_tools import __version__, __mitimroot__ from IPython import embed +from mitim_modules.maestro.utils.EPEDbeat import eped_beat from mitim_modules.maestro.utils.TRANSPbeat import transp_beat from mitim_modules.maestro.utils.PORTALSbeat import portals_beat -from mitim_modules.maestro.utils.EPEDbeat import eped_beat +from mitim_modules.maestro.utils.LENGYELbeat import lengyel_beat from mitim_modules.maestro.utils.MAESTRObeat import creator_from_eped, creator_from_parameterization, creator from mitim_modules.maestro.utils.MAESTRObeat import beat as beat_generic @@ -110,6 +111,9 @@ def define_beat(self, beat, initializer = None, cold_start = False): elif beat == 'eped': print(f'\n- Beat {self.counter_current}: EPED ******************************* {timeBeginning.strftime("%Y-%m-%d %H:%M:%S")}') self.beats[self.counter_current] = eped_beat(self) + elif beat == 'lengyel': + print(f'\n- Beat {self.counter_current}: LENGYEL ******************************* {timeBeginning.strftime("%Y-%m-%d %H:%M:%S")}') + self.beats[self.counter_current] = lengyel_beat(self) # Access current beat easily self.beat = self.beats[self.counter_current] @@ -315,10 +319,14 @@ def _plot_beats(self, fn, num_beats = 2, only_beats = None, full_plot = True): if only_beats is None or only_beats == beat.name: print(f'\t- Plotting beat #{counter}...') - log_file = self.folder_logs / f'plot_{counter}.log' if (not self.terminal_outputs) else None - with LOGtools.conditional_log_to_file(write_log=not ENABLE_EMBED,log_file=log_file): - msg = beat.plot(fn = self.fn, counter = i, full_plot = full_plot) - print(msg) + try: + log_file = self.folder_logs / f'plot_{counter}.log' if (not self.terminal_outputs) else None + with LOGtools.conditional_log_to_file(write_log=not ENABLE_EMBED,log_file=log_file): + msg = beat.plot(fn = self.fn, counter = i, full_plot = full_plot) + print(msg) + except FileNotFoundError: + print(f'\t\t- Could not plot beat #{counter} because some files are missing', typeMsg = 'w') + def _plot_results(self, fn): diff --git a/src/mitim_modules/maestro/scripts/plot_maestro.py b/src/mitim_modules/maestro/scripts/plot_maestro.py index 5c53c1fd..7bf02542 100644 --- a/src/mitim_modules/maestro/scripts/plot_maestro.py +++ b/src/mitim_modules/maestro/scripts/plot_maestro.py @@ -3,6 +3,8 @@ from mitim_tools.misc_tools import GRAPHICStools, IOtools, GUItools from mitim_tools.opt_tools import STRATEGYtools from mitim_tools.misc_tools.utils import remote_tools +from mitim_tools.plasmastate_tools.utils import state_plotting +from mitim_tools.misc_tools.LOGtools import printMsg as print from pathlib import Path from IPython import embed @@ -115,9 +117,11 @@ def main(): colors = GRAPHICStools.listColors() + ms = [] x, scripts = [], [] x0, scripts0 = [], [] + ps_final = [] for i,folder in enumerate(folders): m, ps, ps_lab = MAESTROplot.plotMAESTRO(folder, fn = fn, num_beats=beats, only_beats = only, full_plot = full) ms.append(m) @@ -132,6 +136,16 @@ def main(): if len(x0) > len(x): x = x0 scripts = scripts0 + + ps_final.append(m.final_state) + + # Plot all profiles + if any(ps_final): + figsProfs = state_plotting.add_figures(fn,fnlab_pre = "MAESTRO Profiles ALL - ", tab_color=4) + state_plotting.plotAll(ps_final, figs=figsProfs) + else: + print("No final profiles to plot, no simulation finished", typeMsg="w") + if len(folders) > 1: for let in ['A','B']: diff --git a/src/mitim_modules/maestro/scripts/run_maestro.py b/src/mitim_modules/maestro/scripts/run_maestro.py index a3934e56..95ce9655 100644 --- a/src/mitim_modules/maestro/scripts/run_maestro.py +++ b/src/mitim_modules/maestro/scripts/run_maestro.py @@ -7,6 +7,7 @@ from mitim_modules.maestro.utils import TRANSPbeat, PORTALSbeat from mitim_tools.misc_tools.IOtools import mitim_timer from mitim_tools.misc_tools import PLASMAtools +from mitim_tools.misc_tools.LOGtools import printMsg as print from IPython import embed def profiles_postprocessing_fun(file_profs, lumpImpurities = True, enforce_same_density_gradients = True): @@ -110,7 +111,7 @@ def parse_maestro_nml(file_path): beat_namelists = {} - for beat_type in ["eped","eped_initializer", "transp", "transp_soft", "portals", "portals_soft"]: + for beat_type in ["eped","eped_initializer", "transp", "transp_soft", "portals", "portals_soft", "lengyel"]: if f"{beat_type}_beat" in maestro_namelist["maestro"]: @@ -139,8 +140,6 @@ def parse_maestro_nml(file_path): beat_namelist = maestro_namelist["maestro"][f"{beat_type}_beat"][f"{beat_type}_namelist"] - - # *************************************************************************** # Nothin yet # *************************************************************************** @@ -165,7 +164,8 @@ def parse_maestro_nml(file_path): beat_namelist = maestro_namelist["maestro"]["eped_beat"]["eped_namelist"] else: - raise ValueError(f"[MITIM] {beat_type} beat not found in the MAESTRO namelist") + beat_namelist = None + print(f"[MITIM] {beat_type} beat not found in the MAESTRO namelist", typeMsg='w') beat_namelists[beat_type] = beat_namelist @@ -220,6 +220,8 @@ def run_maestro_local( label_beat = "eped" elif maestro_beats["beats"][0] in ["portals", "portals_soft"]: label_beat = "portals" + elif maestro_beats["beats"][0] in ["lengyel"]: + label_beat = "lengyel" m.define_beat(label_beat, initializer=None if creator_added else parameters_initialize["initializer"]) diff --git a/src/mitim_modules/maestro/utils/EPEDbeat.py b/src/mitim_modules/maestro/utils/EPEDbeat.py index 3504372c..15ae7940 100644 --- a/src/mitim_modules/maestro/utils/EPEDbeat.py +++ b/src/mitim_modules/maestro/utils/EPEDbeat.py @@ -30,12 +30,15 @@ def prepare( BetaN = None, # Force this BetaN (e.g. at creator stage), otherwise from the profiles_current Tesep_keV = None, # Force this Te at the separatrix, otherwise from the profiles_current nesep_20 = None, # Force this ne at the separatrix, otherwise from the profiles_current - corrections_set = {}, # Force these inputs to the NN (e.g. exact delta, Rmajor, etc) + corrections_set = None, # Force these inputs to the NN (e.g. exact delta, Rmajor, etc) ptop_multiplier = 1.0, # Multiplier for the ptop, useful for sensitivity studies TioverTe = 1.0, # Ratio of Ti/Te at the top of the pedestal **kwargs ): + if corrections_set is None: + corrections_set = {} + if nn_location is not None: print(f'\t- Choice of EPED: NN from {IOtools.clipstr(nn_location)}', typeMsg='i') @@ -348,6 +351,7 @@ def _run(self, loopBetaN = 1, minimum_relative_change_in_x=0.005, store_scan = F scan_results[key]['wtop_psipol'] = np.array(scan_results[key]['wtop_psipol']) scan_results[key]['value'] = np.array(scan_results[key]['value']) + self.nn.force_within_range = None # Do not throw warnings during the scan scan_results[key]['ptop_kPa_nominal'], scan_results[key]['wtop_psipol_nominal'] = self.nn(*inputs_to_eped) # --------------------------------- @@ -472,8 +476,8 @@ def plot(self, fn = None, counter = 0, full_plot = True): fig = fn.add_figure(label='EPED', tab_color=counter) axs = fig.subplot_mosaic( """ - ABCDH - AEFGI + ABCDHJ + AEFGIK """ ) axs = [ ax for ax in axs.values() ] diff --git a/src/mitim_modules/maestro/utils/LENGYELbeat.py b/src/mitim_modules/maestro/utils/LENGYELbeat.py new file mode 100644 index 00000000..7353731f --- /dev/null +++ b/src/mitim_modules/maestro/utils/LENGYELbeat.py @@ -0,0 +1,246 @@ +import os +import numpy as np +import copy +from mitim_tools.gacode_tools import PROFILEStools +from mitim_tools.eped_tools import EPEDtools +from mitim_tools.misc_tools import IOtools, GRAPHICStools, GUItools +from mitim_tools.surrogate_tools import NNtools +from mitim_tools.popcon_tools import FunctionalForms +from mitim_tools.misc_tools.LOGtools import printMsg as print +from mitim_modules.maestro.utils.MAESTRObeat import beat +from mitim_modules.powertorch.utils import CALCtools +from mitim_tools.simulation_tools.physics.LENGYELtools import Lengyel +from IPython import embed + +# <> Function to interpolate a curve <> +from mitim_tools.misc_tools.MATHtools import extrapolateCubicSpline as interpolation_function +from sympy import im + +def element_to_lengyel(symbol): + + import periodictable as pt + e = pt.elements.symbol(symbol) # 'W' + + name = e.name # 'tungsten' + charge = e.number # 74 + mass = e.mass # 183.84 + + return name[0].upper() + name[1:].lower(), charge, mass # 'Tungsten', 74, 183.84 + +class lengyel_beat(beat): + + def __init__(self, maestro_instance, folder_name = None): + super().__init__(maestro_instance, beat_name = 'lengyel', folder_name = folder_name) + + def prepare(self, *args, radas_dir = None, seed_impurity_species = None, fixed_impurity_species = None, rhotop=None, **kwargs): + + self.rhotop = rhotop + + if radas_dir is not None: + radas_dir_env = radas_dir + else: + radas_dir_env = os.getenv("RADAS_DIR") + + print('\t- Using provided RADAS_DIR for Lengyel beat preparation:', radas_dir_env) + + # Initialize Lengyel object + self.l = Lengyel() + + # Use seed impurity species from maestro namelist + seed_impurity_symbol = seed_impurity_species["name"] + seed_impurity_ratio_sep_top = seed_impurity_species["ratio_sep_top"] + seed_impurity_name, seed_impurity_Z, seed_impurity_A = element_to_lengyel( seed_impurity_symbol ) + + # High-Z impurity search + fixed_impurity_symbol = fixed_impurity_species + fixed_impurity_name, fixed_impurity_Z, fixed_impurity_A = element_to_lengyel( fixed_impurity_symbol ) + + try: + i_W = np.where(self.profiles_current.profiles['name']==fixed_impurity_symbol)[0][0] + except IndexError: + raise ValueError(f"[MAESTRO][LENGYELbeat] The high-Z impurity species '{fixed_impurity_symbol}' was not found in the input.gacode profiles; please ensure it is present to keep its concentration fixed during the Lengyel beat.") + + fixed_impurity_weights = self.profiles_current.derived['fi_vol'][i_W] + + # Prepare Lengyel with default inputs and changes from GACODE + self.l.prep( + radas_dir = radas_dir_env, + input_gacode = self.profiles_current, + ) + + # ---------------------------------------------------- + # To pass to the run + # ---------------------------------------------------- + + self.lengyel_args = { + 'seed_impurity_species': [ seed_impurity_name ], + 'seed_impurity_weights': [ 1.0 ], + 'fixed_impurity_species': [fixed_impurity_name], + 'fixed_impurity_weights': [fixed_impurity_weights] + } + + # ---------------------------------------------------- + # Other impurity information for post-processing + # ---------------------------------------------------- + self.seed_impurity_enrichment = seed_impurity_ratio_sep_top + + self.seed_impurity_symbol = seed_impurity_symbol + self.seed_impurity_Z = seed_impurity_Z + self.seed_impurity_A = seed_impurity_A + + self.fixed_impurity_symbol = fixed_impurity_symbol + + self._inform() + + def run(self, *args, **kwargs): + + # Run Lengyel standalone + self.l.run( + self.folder, + cold_start=True, # It is so cheap that, if I have come to the run() command, I'll just repeat + **self.lengyel_args + ) + + # Grab important parameters from the inputs + impurity_name = self.lengyel_args['seed_impurity_species'][0] # Assume only one seed impurity in Lengyel run + impurity_symbol = self.seed_impurity_symbol + impurity_Z = self.seed_impurity_Z + impurity_A = self.seed_impurity_A + + # Grab important parameters from the outputs + Tesep = float(self.l.results['separatrix_electron_temp'].split()[0])*1E-3 + + fZ_sep = self.l.results['impurity_fraction']['seed_impurity'][impurity_name] + fZ_top = fZ_sep / self.seed_impurity_enrichment + + # ------------------------------------------------ + # Modify input.gacode + # ------------------------------------------------ + print(f'\t- Applying Lengyel outputs to profiles:') + p = copy.deepcopy(self.profiles_current) + + # Modify temperature profiles + _modify_temperatures(p, Tesep, self.rhotop) + + # Modify impurity density profile + + # Find impurity index: if I have transp beat before Lengyel, I know the order of the impuritites + if "impurity_order_transp" in self.maestro_instance.parameters_trans_beat: + # Impurities ordered as in TRANSP + impurities_in_transp = list(self.maestro_instance.parameters_trans_beat['impurity_order_transp'].keys()) + + # Do not consider the high-Z impurity, which is fixed + impurities_in_transp.remove( self.fixed_impurity_symbol ) + + # Get index of the FIRST (unique for now) seed impurity in TRANSP + i_Z = self.maestro_instance.parameters_trans_beat['impurity_order_transp'][impurities_in_transp[0]] + # Else, assume impurity is in position 3 (after D and He) + else: + i_Z = 3 + print(f"\t\t- No impurity order from TRANSP beat found, assuming impurity '{impurity_symbol}' is in position #{i_Z} in input.gacode", typeMsg='w') + + _modify_impurity_density(p, impurity_symbol, impurity_Z, impurity_A, fZ_sep, fZ_top, self.rhotop, i_Z = i_Z) + + # Enforce quasineutrality + p.enforceQuasineutrality() + + # Check if the plasma just had too much impurity + if p.profiles['ni(10^19/m^3)'].min() < 0: + raise ValueError(f"[MAESTRO][LENGYELbeat] After applying Lengyel outputs, negative ions densities were found in input.gacode.lengyel; please check the impurity concentrations and/or the Lengyel settings.") + + # Write modified input.gacode.lengyel + p.write_state(file=self.folder / 'input.gacode.lengyel') + + def finalize(self, *args, **kwargs): + + self.profiles_output = PROFILEStools.gacode_state(self.folder / 'input.gacode.lengyel') + + self.profiles_output.write_state(file=self.folder_output / 'input.gacode') + + # ----------------------------------------------------------------------------------------------------------------------- + # MAESTRO interface + # ----------------------------------------------------------------------------------------------------------------------- + + def _inform(self, *args, **kwargs): + + # From a previous EPED beat, grab the rhotop + if 'rhotop' in self.maestro_instance.parameters_trans_beat: + self.rhotop = self.maestro_instance.parameters_trans_beat['rhotop'] + print(f"\t\t- Using previous rhotop: {self.rhotop}") + + def _inform_save(self, *args, **kwargs): + + # If I have run Lengyel, I cannot reuse surrogate data #TODO: Maybe not always true? + self.maestro_instance.parameters_trans_beat['portals_surrogate_data_file'] = None + +def _modify_temperatures(p, Tesep, rhotop): + + print(f'\t\t* Setting electron and ion temperature at separatrix to {Tesep*1E3:.1f} eV') + + if rhotop is None: + print('\t\t- No rhotop available at this beat, scaling the entire profile uniformly') + p.profiles['te(keV)'] *= Tesep / p.profiles['te(keV)'][-1] + p.profiles['ti(keV)'] *= Tesep / p.profiles['ti(keV)'][-1, :] + else: + print(f'\t\t- Using rhotop = {rhotop:.3f} to scale temperature profiles only from rhotop to the new separatrix value') + + _scale_quadratic(p, p.profiles['te(keV)'], rhotop, Tesep) + for ion in range(len(p.profiles['ti(keV)'][0, :])): + _scale_quadratic(p, p.profiles['ti(keV)'][:,ion], rhotop, Tesep) + + +def _modify_impurity_density(p, impurity_name, impurity_Z, impurity_A, fZ_sep, fZ_top, rhotop, i_Z): + + print(f'\t\t* Setting impurity "{impurity_name}" (Z={impurity_Z}, A={impurity_A}), at ion position #{i_Z}, density at separatrix to {fZ_top = :.1e}') + + p.profiles['z'][i_Z] = impurity_Z + p.profiles['mass'][i_Z] = impurity_A + p.profiles['name'][i_Z] = impurity_name[:2] + + if rhotop is None: + print('\t\t- No rhotop available at this beat, scaling the entire impurity density profile uniformly by the top (after applying enrichment) value, exact from ne profile') + p.profiles['ni(10^19/m^3)'][:, i_Z] = fZ_top * p.profiles['ne(10^19/m^3)'] + + else: + print(f'\t\t- Using rhotop = {rhotop:.3f} to scale impurity density profiles') + + # First, scale impurity density profile entirely by the desired top value + nZ_top_new = fZ_top * p.profiles['ne(10^19/m^3)'][-1] + ix = np.argmin(np.abs(p.profiles['rho(-)'] - rhotop)) + nZ_top_old = p.profiles['ni(10^19/m^3)'][ix, i_Z] + + p.profiles['ni(10^19/m^3)'][:, i_Z] *= nZ_top_new / nZ_top_old + + # Then modify the edge such that it goes to nZ_sep (Apply quadratic scaling from rhotop to separatrix) + nZ_sep = fZ_sep * p.profiles['ne(10^19/m^3)'][-1] + _scale_quadratic(p, p.profiles['ni(10^19/m^3)'][:, i_Z], rhotop, nZ_sep) + +def _scale_quadratic(p, var, rhotop, val_sep, plotYN=False): + ''' + I use a quadratic scaling from rhotop to separatrix instead of a linear one, to have a smoother transition. + ''' + + var_orig = copy.deepcopy(var) + + ix = np.argmin(np.abs(p.profiles['rho(-)'] - rhotop)) + factor_array = np.ones_like( p.profiles['rho(-)'] ) + + # Create a non-linear scaling that starts,io9slowly and accelerates + n_points = len(p.profiles['rho(-)']) - ix + t = np.linspace(0, 1, n_points) # Normalized parameter from 0 to 1 + factor_array[ix:] = 1.0 + (val_sep / var_orig[-1] - 1.0) * t**2 # Quadratic profile + + var *= factor_array + + if plotYN: + import matplotlib.pyplot as plt + fig, axs= plt.subplots(nrows=2) + ax = axs[0] + ax.plot( p.profiles['rho(-)'], var_orig, 'o-', label='Original' ) + ax.plot( p.profiles['rho(-)'], var, 'o-', label='Modified' ) + ax = axs[1] + ax.plot( p.profiles['rho(-)'], var/var_orig, 'o-', label='Scaling factor' ) + plt.show() + + embed() + \ No newline at end of file diff --git a/src/mitim_modules/maestro/utils/MAESTROplot.py b/src/mitim_modules/maestro/utils/MAESTROplot.py index 685c639b..dc8c4f30 100644 --- a/src/mitim_modules/maestro/utils/MAESTROplot.py +++ b/src/mitim_modules/maestro/utils/MAESTROplot.py @@ -13,6 +13,7 @@ from mitim_modules.maestro.utils.TRANSPbeat import transp_beat from mitim_modules.maestro.utils.PORTALSbeat import portals_beat from mitim_modules.maestro.utils.EPEDbeat import eped_beat +from mitim_modules.maestro.utils.LENGYELbeat import lengyel_beat MARKERSIZE = 1 LW = 1.0 @@ -21,7 +22,7 @@ def grabMAESTRO(folder): # Find beat results from folders folder_beats = Path(folder) / 'Beats' - beats = sorted([item.name for item in folder_beats.glob('*') if not item.name.startswith(".")], key=lambda x: int(x.split('_')[1])) + beats = sorted([item.name for item in folder_beats.glob('*') if item.name.startswith("Beat")], key=lambda x: int(x.split('_')[1])) beat_types = [] for beat in range(len(beats)): @@ -31,6 +32,8 @@ def grabMAESTRO(folder): beat_types.append('portals') elif (folder_beats / f'{beats[beat]}' / 'run_eped').exists(): beat_types.append('eped') + elif (folder_beats / f'{beats[beat]}' / 'run_lengyel').exists(): + beat_types.append('lengyel') # First initializer beat_initializer = None @@ -47,6 +50,13 @@ def grabMAESTRO(folder): for i,beat in enumerate(beat_types): m.define_beat(beat, initializer = beat_initializer if i == 0 else None) + # Add final if exists + folder_output = Path(folder) / 'Outputs' + if (folder_output / 'input.gacode_final').exists(): + m.final_state = PROFILEStools.gacode_state(folder_output / 'input.gacode_final') + else: + m.final_state = None + return m def plotMAESTRO(folder, fn = None, num_beats = 2, only_beats = None, full_plot = True): @@ -90,6 +100,8 @@ def plot_results(self, fn): key = f'PORTALS b#{i+1}' elif isinstance(beat, eped_beat): key = f'EPED b#{i+1}' + elif isinstance(beat, lengyel_beat): + key = f'Lengyel b#{i+1}' objs[key] = profs @@ -123,8 +135,8 @@ def plot_results(self, fn): fig = fn.add_figure(label='MAESTRO init', tab_color=2) axs = fig.subplot_mosaic( """ - ABCDH - AEFGI + ABCDHK + AEFGIJ """ ) axs = [ ax for ax in axs.values() ] @@ -153,8 +165,8 @@ def plot_results(self, fn): fig = fn.add_figure(label=f'{label} {i}->{i+1}', tab_color=2) axs = fig.subplot_mosaic( """ - ABCDH - AEFGI + ABCDHJ + AEFGIK """ ) axs = [ ax for ax in axs.values() ] @@ -171,8 +183,8 @@ def plot_results(self, fn): fig = fn.add_figure(label=f'{label} {0}->{len(keys)}', tab_color=2) axs = fig.subplot_mosaic( """ - ABCDH - AEFGI + ABCDHJ + AEFGIK """ ) axs = [ ax for ax in axs.values() ] @@ -233,7 +245,7 @@ def plot_special_quantities(ps, ps_lab, axs, color='b', label = '', legYN=True): nu_ne.append(p.derived['ne_peaking0.2']) q95.append(p.derived['q95']) q0.append(p.derived['q0']) - xsaw.append(p.derived['rho_saw']) + xsaw.append(p.derived['roa_saw']) p90.append(np.interp(0.9,p.profiles['rho(-)'],p.derived['pthr_manual'])) def _special(ax,x): @@ -343,7 +355,7 @@ def _special(ax,x): ax.axhline(y=1, color = 'k', lw = 2, ls = '--') if legYN: ax.legend() - ax.set_ylim(bottom = 0) + #ax.set_ylim(bottom = 0) ax.set_xticklabels([]) @@ -351,7 +363,7 @@ def _special(ax,x): ax = axs['J'] ax.plot(x, xsaw, '-s', color=color, markersize=7, lw = 1) - ax.set_ylabel('Inversion radius (rho)') + ax.set_ylabel('Inversion radius (roa)') GRAPHICStools.addDenseAxis(ax) ax.set_ylim([0,1]) diff --git a/src/mitim_modules/maestro/utils/PORTALSbeat.py b/src/mitim_modules/maestro/utils/PORTALSbeat.py index b3c43074..fe815042 100644 --- a/src/mitim_modules/maestro/utils/PORTALSbeat.py +++ b/src/mitim_modules/maestro/utils/PORTALSbeat.py @@ -106,7 +106,11 @@ def run(self, **kwargs): self.mitim_bo = STRATEGYtools.MITIM_BO(portals_fun, seed = self.maestro_instance.master_seed, cold_start = cold_start, askQuestions = False) - if self.use_previous_surrogate_data and self.try_flux_match_only_for_first_point and self.folder_starting_point is not None: + if self.use_previous_surrogate_data and \ + self.try_flux_match_only_for_first_point and \ + self.folder_starting_point is not None and \ + ('portals_surrogate_data_file' in self.maestro_instance.parameters_trans_beat) and \ + self.maestro_instance.parameters_trans_beat['portals_surrogate_data_file'] is not None: # PORTALS just with one point portals_fun.optimization_options['initialization_options']['initial_training'] = 1 diff --git a/src/mitim_modules/maestro/utils/TRANSPbeat.py b/src/mitim_modules/maestro/utils/TRANSPbeat.py index 14936211..6fb965e3 100644 --- a/src/mitim_modules/maestro/utils/TRANSPbeat.py +++ b/src/mitim_modules/maestro/utils/TRANSPbeat.py @@ -1,6 +1,7 @@ import os import shutil import copy +from typing import OrderedDict from mitim_tools.transp_tools import CDFtools from mitim_tools.misc_tools import IOtools from mitim_tools.gacode_tools import PROFILEStools @@ -64,7 +65,6 @@ def prepare( self.time_end = self.time_diffusion + flattop_window # End self.timeAC = self.time_end - 0.001 if extractAC else None # Time to extract TORIC and NUBEAM files - # Write TRANSP from profiles times = [self.time_transition,self.time_end+1.0] self.transp = self.profiles_current.to_transp( @@ -147,7 +147,7 @@ def run(self, **kwargs): def finalize(self, force_auxiliary_heating_at_output = {'Pe': None, 'Pi': None}, **kwargs): # Copy to outputs - try: + try: shutil.copy2(self.folder / f"{self.shot}{self.runid}TR.DAT", self.folder_output) shutil.copy2(self.folder / f"{self.shot}{self.runid}.CDF", self.folder_output) shutil.copy2(self.folder / f"{self.shot}{self.runid}tr.log", self.folder_output) @@ -168,6 +168,17 @@ def finalize(self, force_auxiliary_heating_at_output = {'Pe': None, 'Pi': None}, shutil.copy2(self.folder / f"{cdf_prefix}.CDF", self.folder_output / f"{self.shot}{self.runid}.CDF") shutil.copy2(self.folder / f"{cdf_prefix}tr.log", self.folder_output / f"{self.shot}{self.runid}tr.log") + # Remove any existing files in the output folder (to avoid multiple CDFs) + for cdf_file in self.folder_output.glob("*.CDF"): + if cdf_file.name != f"{self.shot}{self.runid}.CDF": + os.remove(cdf_file) + for trlog_file in self.folder_output.glob("*tr.log"): + if trlog_file.name != f"{self.shot}{self.runid}tr.log": + os.remove(trlog_file) + for trdat_file in self.folder_output.glob("*TR.DAT"): + if trdat_file.name != f"{self.shot}{self.runid}TR.DAT": + os.remove(trdat_file) + # Extract output cdf_results = CDFtools.transp_output(self.folder_output / f"{self.shot}{self.runid}.CDF") @@ -299,7 +310,25 @@ def _additional_operations_add_initialization(self, machine_initialization = 'CM # ----------------------------------------------------------------------------------------------------------------------- def _inform_save(self, *args, **kwargs): - self.maestro_instance.parameters_trans_beat['portals_surrogate_data_file'] = None # If I have run TRANSP, I cannot reuse surrogate data #TODO: Maybe not always true? + c, _ = self.grab_output() + + # Grab the oder of user-specified impuritites in the TRANSP ions list + + transp_impurities = c.nZs.keys() + profiles_species = [i['N'] for i in self.profiles_output.Species] + + impurity_order_transp = OrderedDict() + for z in transp_impurities: + for i,spec in enumerate(profiles_species): + if spec == z: + impurity_order_transp[spec] = i + break + + self.maestro_instance.parameters_trans_beat['impurity_order_transp'] = impurity_order_transp + + # If I have run TRANSP, I cannot reuse surrogate data #TODO: Maybe not always true? + + self.maestro_instance.parameters_trans_beat['portals_surrogate_data_file'] = None # ----------------------------------------------------------------------------------------------------------------------- # Defaults to help MAESTRO diff --git a/src/mitim_tools/gs_tools/GEQtools.py b/src/mitim_tools/gs_tools/GEQtools.py index d6a349ef..5595344e 100644 --- a/src/mitim_tools/gs_tools/GEQtools.py +++ b/src/mitim_tools/gs_tools/GEQtools.py @@ -27,7 +27,10 @@ class MITIMgeqdsk: def __init__(self, filename): self.g = megpy.Equilibrium() - self.g.read_geqdsk(f_path=filename) + try: + self.g.read_geqdsk(f_path=filename) + except ValueError: + raise ValueError("-> MITIMgeqdsk: Problem reading g-eqdsk file ", filename) self.g.add_derived(incl_fluxsurfaces=True, analytic_shape=True, incl_B=True) # Extra derivations in MITIM diff --git a/src/mitim_tools/misc_tools/IOtools.py b/src/mitim_tools/misc_tools/IOtools.py index c431b689..ec33db9a 100644 --- a/src/mitim_tools/misc_tools/IOtools.py +++ b/src/mitim_tools/misc_tools/IOtools.py @@ -945,7 +945,7 @@ def findFileByExtension( if len(allfiles) > 1: # print(allfiles) if not ForceFirst: - raise Exception("More than one file with same extension in the folder!") + raise Exception("More than one file with same extension in the folder: ", fpath) else: allfiles = [allfiles[0]] diff --git a/src/mitim_tools/plasmastate_tools/MITIMstate.py b/src/mitim_tools/plasmastate_tools/MITIMstate.py index 81eba509..a28a3403 100644 --- a/src/mitim_tools/plasmastate_tools/MITIMstate.py +++ b/src/mitim_tools/plasmastate_tools/MITIMstate.py @@ -355,8 +355,10 @@ def derive_quantities_full(self, mi_ref=None, rederiveGeometry=True): if self.profiles["q(-)"].min() > 1.0: self.derived["rho_saw"] = np.nan + self.derived["roa_saw"] = np.nan else: self.derived["rho_saw"] = np.interp(1.0, self.profiles["q(-)"], self.profiles["rho(-)"]) + self.derived["roa_saw"] = np.interp(1.0, self.profiles["q(-)"], self.derived["roa"]) # --------- Geometry (only if it doesn't exist or if I ask to recalculate) @@ -2048,11 +2050,11 @@ def plotRelevant(self, axs = None, color = 'b', label ='', lw = 1, ms = 1): fig = plt.figure() axs = fig.subplot_mosaic( """ - ABCDH - AEFGI + ABCDHJ + AEFGIK """ ) - axs = [axs['A'], axs['B'], axs['C'], axs['D'], axs['E'], axs['F'], axs['G'], axs['H'], axs['I']] + axs = [axs['A'], axs['B'], axs['C'], axs['D'], axs['E'], axs['F'], axs['G'], axs['H'], axs['I'], axs['J'], axs['K']] # ---------------------------------- # Equilibria @@ -2138,38 +2140,26 @@ def plotRelevant(self, axs = None, color = 'b', label ='', lw = 1, ms = 1): # Powers # ---------------------------------- - # RF + # RF and Ohmic ax = axs[5] - ax.plot(self.profiles['rho(-)'], self.profiles['qrfe(MW/m^3)'], '-o', markersize=ms, lw = lw, label=label+', e', color=color) - ax.plot(self.profiles['rho(-)'], self.profiles['qrfi(MW/m^3)'], '--*', markersize=ms, lw = lw, label=label+', i', color=color) + ax.plot(self.profiles['rho(-)'], self.profiles['qrfe(MW/m^3)'], '-o', markersize=ms, lw = lw, label=label+', ICH e', color=color) + ax.plot(self.profiles['rho(-)'], self.profiles['qrfi(MW/m^3)'], '--*', markersize=ms, lw = lw, label=label+', ICH i', color=color) + ax.plot(self.profiles['rho(-)'], self.profiles['qohme(MW/m^3)'], '-.s', markersize=ms, lw = lw, label=label+', Ohmic', color=color) ax.set_xlabel("$\\rho_N$") - ax.set_ylabel("$P_{ich}$ (MW/m$^3$)") + ax.set_ylabel("$P$ (MW/m$^3$)") #ax.set_ylim(bottom = 0) ax.set_xlim(0,1) ax.legend(prop={'size':8}) GRAPHICStools.addDenseAxis(ax) - ax.set_title("ICH Power Deposition") - - # Ohmic - ax = axs[6] - - ax.plot(self.profiles['rho(-)'], self.profiles['qohme(MW/m^3)'], '-o', markersize=ms, lw = lw, label=label, color=color) - - ax.set_xlabel("$\\rho_N$") - ax.set_ylabel("$P_{oh}$ (MW/m$^3$)") - #ax.set_ylim(bottom = 0) - ax.set_xlim(0,1) - ax.legend(prop={'size':8}) - GRAPHICStools.addDenseAxis(ax) - ax.set_title("Ohmic Power Deposition") + ax.set_title("Power Deposition") # ---------------------------------- # Heat fluxes # ---------------------------------- - ax = axs[7] + ax = axs[6] ax.plot(self.profiles['rho(-)'], self.derived['qe_MWm2'], '-o', markersize=ms, lw = lw, label=label+', e', color=color) ax.plot(self.profiles['rho(-)'], self.derived['qi_MWm2'], '--*', markersize=ms, lw = lw, label=label+', i', color=color) @@ -2186,7 +2176,7 @@ def plotRelevant(self, axs = None, color = 'b', label ='', lw = 1, ms = 1): # Dynamic targets # ---------------------------------- - ax = axs[8] + ax = axs[7] ax.plot(self.profiles['rho(-)'], self.derived['qrad'], '-o', markersize=ms, lw = lw, label=label+', rad', color=color) ax.plot(self.profiles['rho(-)'], self.profiles['qei(MW/m^3)'], '--*', markersize=ms, lw = lw, label=label+', exc', color=color) @@ -2201,6 +2191,77 @@ def plotRelevant(self, axs = None, color = 'b', label ='', lw = 1, ms = 1): GRAPHICStools.addDenseAxis(ax) ax.set_title("Dynamic Targets") + # ---------------------------------- + # Ions: Main + # ---------------------------------- + + ls = GRAPHICStools.listmarkersLS() + + ax = axs[8] + cont = 0 + addedSpecies = [] + sumFi = np.zeros(self.profiles['rho(-)'].shape) + if self.Dion is not None: + ax.plot(self.profiles['rho(-)'], self.derived['fi'][:,self.Dion], ls[cont], markersize=ms, lw = lw, label=label+', D', color=color) + addedSpecies.append(self.Dion) + sumFi += self.derived['fi'][:,self.Dion] + cont += 1 + if self.Tion is not None: + ax.plot(self.profiles['rho(-)'], self.derived['fi'][:,self.Tion], ls[cont], markersize=ms, lw = lw, label=label+', T', color=color) + addedSpecies.append(self.Tion) + sumFi += self.derived['fi'][:,self.Tion] + cont += 1 + if 'Mion' in self.__dict__ and self.Mion is not None: + ax.plot(self.profiles['rho(-)'], self.derived['fi'][:,self.Mion], ls[cont], markersize=ms, lw = lw, label=label+', Main', color=color) + addedSpecies.append(self.Mion) + sumFi += self.derived['fi'][:,self.Mion] + cont += 1 + + ax.plot(self.profiles['rho(-)'], sumFi, ls[cont], markersize=ms, lw = lw*2.0, label=label+', Sum Main', color=color) + + ax.set_xlabel("$\\rho_N$") + ax.set_ylabel("Concentration ($f_i$)") + ax.set_ylim([0,1]) + ax.set_xlim(0,1) + ax.legend(prop={'size':8}) + GRAPHICStools.addDenseAxis(ax) + ax.set_title("Main Ions") + + # ---------------------------------- + # Ions: Others + # ---------------------------------- + + ax = axs[9] + cont = 0 + for i in range(len(self.Species)): + if i not in addedSpecies: + ax.plot(self.profiles['rho(-)'], self.derived['fi'][:,i], ls[cont], markersize=ms, lw = lw, label=label+f', {self.profiles["name"][i]}', color=color) + cont += 1 + + ax.set_xlabel("$\\rho_N$") + ax.set_ylabel("Concentration ($f_i$)") + #ax.set_ylim(bottom = 0) + ax.set_yscale('log') + ax.set_xlim(0,1) + ax.legend(prop={'size':8}) + GRAPHICStools.addDenseAxis(ax) + ax.set_title("Impurity Ions") + + # ---------------------------------- + # Ions: Zeff + # ---------------------------------- + + ax = axs[10] + ax.plot(self.profiles['rho(-)'], self.derived['Zeff'], '-o', markersize=ms, lw = lw, label=label, color=color) + ax.set_xlabel("$\\rho_N$") + ax.set_ylabel("Zeff") + #ax.set_ylim(bottom = 0) + ax.set_xlim(0,1) + ax.legend(prop={'size':8}) + GRAPHICStools.addDenseAxis(ax) + ax.set_title("Zeff Profile") + + def csv(self, file="input.gacode.xlsx"): dictExcel = IOtools.OrderedDict() diff --git a/src/mitim_tools/plasmastate_tools/utils/state_plotting.py b/src/mitim_tools/plasmastate_tools/utils/state_plotting.py index 9559e42f..5e547741 100644 --- a/src/mitim_tools/plasmastate_tools/utils/state_plotting.py +++ b/src/mitim_tools/plasmastate_tools/utils/state_plotting.py @@ -19,7 +19,6 @@ def add_figures(fn, fnlab='', fnlab_pre='', tab_color=None): return figs - def add_axes(figs): fig1, fig2, fig3, fig4, fig5, fig6, fig7 = figs @@ -1337,6 +1336,8 @@ def plotAll(profiles_list, figs=None, extralabs=None, lastRhoGradients=0.89): else: extralab = f"{extralabs[i]}, " + if profiles is None: + continue profiles.plot( axs1=axsProf_1,axs2=axsProf_2,axs3=axsProf_3,axs4=axsProf_4,axsFlows=axsFlows,axs6=axsProf_6,axsImps=axsImps, color=colors[i],legYN=True,extralab=extralab,lsFlows=ls[i],legFlows=i == 0,showtexts=False,lastRhoGradients=lastRhoGradients, diff --git a/src/mitim_tools/simulation_tools/physics/LENGYELtools.py b/src/mitim_tools/simulation_tools/physics/LENGYELtools.py new file mode 100644 index 00000000..bc73fffc --- /dev/null +++ b/src/mitim_tools/simulation_tools/physics/LENGYELtools.py @@ -0,0 +1,140 @@ +import os +from pathlib import Path +import matplotlib.pyplot as plt +from mitim_tools.misc_tools.LOGtools import printMsg as print +from mitim_tools.misc_tools import IOtools, GRAPHICStools +from mitim_tools.gacode_tools import PROFILEStools +from mitim_tools import __mitimroot__ +from IPython import embed + +class Lengyel(): + def __init__( + self + ): + + self.nml_default = Path(__mitimroot__ / 'templates' / 'input.lengyel.controls.yaml') + + # Optional preparation step + def prep( + self, + radas_dir, + input_gacode = None + ): + + # Read default namelist + self.nml = IOtools.read_mitim_yaml(self.nml_default) + + # Change RADAS directory + radas_dir = Path(radas_dir) + if radas_dir.exists(): + self.nml['input']['radas_dir'] = f"PATH:{radas_dir.resolve()}" + else: + raise FileNotFoundError(f"[MITIM] The provided RADAS_DIR '{radas_dir}' does not exist; please do 'radas -c radas_config.yml -s tungsten' with the proper config and impurities") + + # Potentially change some parameters from the input.gacode + + params = { + 'major_radius': ['profiles', 'rcentr(m)', 'm', 0], + 'minor_radius': ['derived', 'a', 'm', None], + 'elongation_psi95': ['derived', 'kappa95', ' ', None], + 'triangularity_psi95': ['derived', 'delta95', ' ', None], + 'magnetic_field_on_axis': ['profiles', 'bcentr(T)', 'T', 0], + 'plasma_current': ['profiles', 'current(MA)', 'MA', 0], + 'ion_mass': ['derived', 'mbg_main', 'amu', None], + 'power_crossing_separatrix': ['derived', 'Psol', 'MW', None], + 'separatrix_electron_density': ['profiles', 'ne(10^19/m^3)', 'e19/m^3', -1], + } + + if input_gacode is not None: + print(f"\t- Populating Lengyel input from provided GACODE profile:") + if isinstance(input_gacode, PROFILEStools.gacode_state): + p = input_gacode + else: + p = PROFILEStools.gacode_state(input_gacode) + + for par in params: + val = p.__dict__[params[par][0]][params[par][1]] + if params[par][3] is not None: + val = val[params[par][3]] + print(f"\t\t* Setting '{par}' to MITIMstate value '{params[par][1]} = {val}'") + self.nml['input'][par] = f'{val}{params[par][2]}' + + def run( + self, + folder, + cold_start = False, + **input_dict + ): + + folder = Path(folder) + if not folder.exists(): + folder.mkdir(parents=True, exist_ok=True) + elif cold_start: + print(f"\t- Cold starting Lengyel run; cleaning folder '{folder}'") + os.system(f'rm -rf {folder}/*') + + # Potentially modify namelist with input_dict + for key in input_dict: + print(f"\t- Setting Lengyel input parameter '{key}' to '{input_dict[key]}'") + self.nml['input'][key] = input_dict[key] + + # Write modified namelist to folder + nml_file = folder / 'input.lengyel.controls.yml' + IOtools.write_mitim_yaml( self.nml, nml_file ) + + # Run + output_file = folder / 'output.lengyel.results.yml' + from extended_lengyel.cli import run_extended_lengyel + run_extended_lengyel( + config_file = nml_file, + output_file = output_file, + ) + + # Read output + self.results = IOtools.read_mitim_yaml( output_file ) + + def run_scan( + self, + folder, + scan_name, + scan_values, + cold_start = False, + plotYN = True, + **input_dict + ): + + folder = Path(folder) + if not folder.exists(): + folder.mkdir(parents=True, exist_ok=True) + + self.results_scan = {} + for val in scan_values: + print(f"\t- Running Lengyel scan '{scan_name}' with value '{val}'") + scan_folder = folder / f"{scan_name}_{val}" + scan_input = input_dict.copy() + scan_input[scan_name] = val + self.run( + folder = scan_folder, + cold_start = cold_start, + **scan_input + ) + self.results_scan[val] = self.results + + # Plot + if plotYN: + plt.ion(); fig, ax = plt.subplots() + for val in scan_values: + res = self.results_scan[val] + ax.plot( + val, + float(res['separatrix_electron_temp'].split()[0]), + 'o', markersize=15 + ) + + ax.set_xlabel(f"{scan_name}") + ax.set_ylabel("Separatrix electron temperature [eV]") + GRAPHICStools.addDenseAxis(ax) + + + + \ No newline at end of file diff --git a/templates/input.lengyel.controls.yaml b/templates/input.lengyel.controls.yaml new file mode 100644 index 00000000..8084cdc9 --- /dev/null +++ b/templates/input.lengyel.controls.yaml @@ -0,0 +1,74 @@ +algorithms: + - calc_magnetic_field_and_safety_factor + - calc_fieldline_pitch_at_omp + - read_atomic_data + - calc_kappa_e0 + - build_CzLINT_for_seed_impurities + - build_mean_charge_for_seed_impurities + - build_CzLINT_for_fixed_impurities + - build_mean_charge_for_fixed_impurities + - calc_momentum_loss_from_cc_fit + - calc_power_loss_from_cc_fit + - calc_electron_temp_from_cc_fit + - run_extended_lengyel_model_with_S_Zeff_and_alphat_correction + - calc_sound_speed_at_target + - calc_target_density + - calc_flux_density_to_pascals_factor + - calc_parallel_to_perp_factor + - calc_ion_flux_to_target + - calc_divertor_neutral_pressure + - calc_heat_flux_perp_to_target +input: + # Set the path to the atomic data directory made by radas. + radas_dir: "PATH:./radas_dir" + # Length along a magnetic fieldline from divertor target to + # divertor entrance. + divertor_parallel_length: 10.0 m + # Length along a magnetic fieldline from divertor target to + # outboard midplane or stagnation point. + parallel_connection_length: 30 m + # Major radius + major_radius: 1.85 m + # Minor radius + minor_radius: 0.57 m + # Elongation at psiN = 0.95 + elongation_psi95: 1.7 + # Triangularity at psiN = 0.95 + triangularity_psi95: 0.3 + # Magnetic field on axis + magnetic_field_on_axis: 12.2 T + # Plasma current + plasma_current: 8.7 MA + # Ratio of Bpol_OMP / Bpol_avg + ratio_of_upstream_to_average_poloidal_field: 4/3 + # Main ion mass (or average mass, i.e. DT = 2.5amu) + ion_mass: 2.5 amu + # List of impurities used for seeding, and their relative concentrations. + seed_impurity_species: ["Neon"] + seed_impurity_weights: [1.0] + # List of background impurities, and their concentration relative to the electron density. + fixed_impurity_species: ["Tungsten"] + fixed_impurity_weights: [1.5e-5] + # Take any nearest ne or ne_tau value for the atomic data, regardless of how close it is + # to the reference value. + rtol_nearest_for_atomic_data: 1e-1 + # Impurity residence time multiplied by electron density + reference_ne_tau: 0.5 ms n20 + # gamma_sh + sheath_heat_transmission_factor: 7.5 + # Angle of incidence between magnetic fieldline and divertor target + target_angle_of_incidence: 3 degree + # Ratio between lambda_INT and lambda_q + divertor_broadening_factor: 3.0 + # Power crossing separatrix + power_crossing_separatrix: 20.9MW + # Fraction of power directed to outer divertor + fraction_of_P_SOL_to_divertor: 2/3 + # Upstream electron density + separatrix_electron_density: 15.0e19/m^3 + # Electron temperature at target + target_electron_temp: 2.34eV + # Set number of iterations for Lengyel model. Should set this large, so that the result doesn't + # change if you keep increasing it. + iterations_for_alphat: 50 + iterations_for_Lengyel_model: 5 \ No newline at end of file diff --git a/templates/namelist.maestro.yaml b/templates/namelist.maestro.yaml index 52c8d7fc..5da8eec9 100644 --- a/templates/namelist.maestro.yaml +++ b/templates/namelist.maestro.yaml @@ -91,13 +91,16 @@ assumptions: maestro: # Sequence of beats - beats: ["transp_soft", "transp", "eped", "portals", "eped", "portals"] + beats: ["transp_soft", "eped", "portals", "eped", "portals"] # Remove intermediate files to avoid heavy MAESTRO simulation folders keep_all_files: true # --------------------------------------------------------------------------- # Each individual beat parameters + # use_default: true/false + # ****_namelist: Parameters that are passed to the run() command + # #TODO # --------------------------------------------------------------------------- eped_beat: @@ -121,7 +124,9 @@ maestro: TioverTe: 1.0 # Ti/Te assumption at the pedestal top eped_initializer_beat: + use_default: false + eped_initializer_namelist: nn_location: $MFEIM_PATH/private_code_mitim/NN_DATA/EPED-NN-SPARC/EPED-NN-MODEL-SPARC.keras norm_location: $MFEIM_PATH/private_code_mitim/NN_DATA/EPED-NN-SPARC/EPED-NN-NORMALIZATION-SPARC.txt @@ -129,9 +134,28 @@ maestro: Bt: 12.2 R: 1.85 a: 0.57 + ptop_multiplier: 1.0 TioverTe: 1.0 + lengyel_beat: + + use_default: false + + lengyel_namelist: + + # If null, it will use the default namelist at: __mitimroot__ / "templates" / "input.lengyel.controls.yaml" + lengyel_namelist_location: null + + # Path to radas directory. If null, it will use $RADAS_DIR environement variable + radas_dir: null + + seed_impurity_species: + name: "Ne" + ratio_sep_top: 1.0 + + fixed_impurity_species: 'W' # It will try find this element in the ions available in input.gacode, keeping their concentration fixed + portals_beat: use_default: false @@ -145,13 +169,11 @@ maestro: # PORTALS parameters (only those parameters to change the main PORTALS namelist) # ------------------------------------------------------------------------------------------------ portals_parameters: + + # Same structure as the PORTALS namelist, here indicate only the parameters to change from default portals_namelist_location solution: predicted_roa: [0.35, 0.45, 0.55, 0.65, 0.75, 0.875, 0.9] keep_full_model_folder: false - exploration_ranges: - ymax: 4.0 - ymin: 1.0 - yminymax_atleast: [null, 4] transport: options: tglf: @@ -172,7 +194,6 @@ maestro: acquisition_options: optimizers: ["sr"] - # Operations to do to mitim_state before passing it to PORTALS initialization_parameters: thermalize_fast: true diff --git a/tests/LENGYEL_workflow.py b/tests/LENGYEL_workflow.py new file mode 100644 index 00000000..7023ba45 --- /dev/null +++ b/tests/LENGYEL_workflow.py @@ -0,0 +1,42 @@ +import os +from mitim_tools.simulation_tools.physics.LENGYELtools import Lengyel +from mitim_tools import __mitimroot__ + +cold_start = True + +(__mitimroot__ / 'tests' / 'scratch').mkdir(parents=True, exist_ok=True) + +folder = __mitimroot__ / "tests" / "scratch" / "lengyel_test" +input_gacode = __mitimroot__ / "tests" / "data" / "input.gacode" + +# Check if RADAS_DIR is set +radas_dir_env = os.getenv("RADAS_DIR") +if radas_dir_env is None: + raise EnvironmentError("[MITIM] The RADAS_DIR environment variable is not set") + +if cold_start and folder.exists(): + os.system(f"rm -r {folder.resolve()}") + +# Initialize Lengyel object +l = Lengyel() + +# Prepare Lengyel with default inputs and changes from GACODE +l.prep( + radas_dir = radas_dir_env, + input_gacode = input_gacode + ) + +# Run Lengyel standalone +l.run( + folder / 'tmp_run', + cold_start=cold_start + ) + +# Run Lengyel scans +l.run_scan( + folder = folder / 'tmp_scan', + scan_name = 'power_crossing_separatrix', + scan_values = ['10MW','20MW','30MW'], + plasma_current = '2.0 MA', # Change any parameter from the default + plotYN = True # Plot results after the scan + )