freegsnke.nonlinear_solve module

Implements the core non-linear solver for the evolutive GS problem. Also handles the linearised evolution capabilites.

Copyright 2025 UKAEA, UKRI-STFC, and The Authors, as per the COPYRIGHT and README files.

This file is part of FreeGSNKE.

FreeGSNKE is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

FreeGSNKE is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

You should have received a copy of the GNU Lesser General Public License along with FreeGSNKE. If not, see <http://www.gnu.org/licenses/>.

class freegsnke.nonlinear_solve.nl_solver(profiles, eq, custom_coil_resist=None, custom_self_ind=None, full_timestep=0.0001, max_internal_timestep=0.0001, automatic_timestep=False, plasma_resistivity=1e-06, plasma_norm_factor=1000.0, blend_hatJ=0, max_mode_frequency=100.0, fix_n_vessel_modes=-1, threshold_dIy_dI=0.2, min_dIy_dI=0.1, mode_removal=True, linearize=True, dIydI=None, target_relative_tolerance_linearization=1e-08, target_dIy=0.001, force_core_mask_linearization=False, verbose=False)

Bases: object

Handles all time-evolution capabilites. Includes interface to use both: - stepper of the linearised problem - stepper for the full non-linear problem

F_function_curr(trial_currents, active_voltage_vec)

Full non-linear system of circuit eqs written as root problem in the vector of current values at time t+dt. Note that the plasma flux at time t+dt is taken to be self.trial_plasma_psi. Iteration consists of: [trial_currents, plasma_flux] -> hatIy1, through calculating plasma distribution hatIy1 -> iterated_currents, through ‘simplified’ circuit eqs Residual: iterated_currents - trial_currents Residual is zero if the pair [trial_currents, plasma_flux] solve the full non-linear problem.

Parameters:

• trial_currents (np.array) – Vector of current values. Same format as self.currents_vec.

• active_voltage_vec (np.array) – Vector of active voltages for the active coils, applied between t and t+dt.

Returns:

Residual in current values. Same format as self.currents_vec.

Return type:

np.array

F_function_curr_GS(trial_currents, active_voltage_vec, rtol_NK)

Full non-linear system of circuit eqs written as root problem in the vector of current values at time t+dt. Note that, differently from self.F_function_curr, here the plasma flux is not imposed, but self-consistently solved for based on the input trial_currents. Iteration consists of: trial_currents -> plasma flux, through static GS [trial_currents, plasma_flux] -> hatIy1, through calculating plasma distribution hatIy1 -> iterated_currents, through ‘simplified’ circuit eqs Residual: iterated_currents - trial_currents Residual is zero if trial_currents solve the full non-linear problem.

Parameters:

• trial_currents (np.array) – Vector of current values. Same format as self.currents_vec.

• active_voltage_vec (np.array) – Vector of active voltages for the active coils, applied between t and t+dt.

• rtol_NK (float) – Relative tolerance to be used in the static GS problem.

Returns:

Residual in current values. Same format as self.currents_vec.

Return type:

np.array

F_function_psi(trial_plasma_psi, active_voltage_vec, rtol_NK)

Full non-linear system of circuit eqs written as root problem in the plasma flux. Note that the flux associated to the metal currents is sourced from self.tokamak_psi. Iteration consists of: [trial_plasma_psi, tokamak_psi] -> hatIy1, by calculating Jtor hatIy1 -> currents(t+dt), through ‘simplified’ circuit eq currents(t+dt) -> iterated_plasma_flux, through static GS Residual: iterated_plasma_flux - trial_plasma_psi Residual is zero if the pair [trial_plasma_psi, tokamak_psi] solve the full non-linear problem.

Parameters:

• trial_plasma_psi (np.array) – Plasma flux values in 1d vector covering full domain of size eq.nx*eq.ny.

• active_voltage_vec (np.array) – Vector of active voltages for the active coils, applied between t and t+dt.

• rtol_NK (float) – Relative tolerance to be used in the static GS problem.

Returns:

Residual in plasma flux, 1d.

Return type:

np.array

__init__(profiles, eq, custom_coil_resist=None, custom_self_ind=None, full_timestep=0.0001, max_internal_timestep=0.0001, automatic_timestep=False, plasma_resistivity=1e-06, plasma_norm_factor=1000.0, blend_hatJ=0, max_mode_frequency=100.0, fix_n_vessel_modes=-1, threshold_dIy_dI=0.2, min_dIy_dI=0.1, mode_removal=True, linearize=True, dIydI=None, target_relative_tolerance_linearization=1e-08, target_dIy=0.001, force_core_mask_linearization=False, verbose=False)

Initializes the time-evolution Object.

Parameters:

• profiles (FreeGSNKE profiles Object) – profiles function of the initial equilibrium. This will be used to set up the linearization used by the linear evolutive solver. It can be changed later by initializing a new set of initial conditions. Note however that, to change either the machine or limiter properties it will be necessary to instantiate a new nl_solver object.

• eq (FreeGSNKE equilibrium Object) – Initial equilibrium. This is used to set the domain/grid properties as well as the machine properties. Furthermore, eq will be used to set up the linearization used by the linear evolutive solver. It can be changed later by initializing a new set of initial conditions. Note however that, to change either the machine or limiter properties it will be necessary to instantiate a new nl_solver object.

• max_mode_frequency (float) – Threshold value used to include/exclude vessel normal modes. Only modes with smaller characteristic frequencies (larger timescales) are retained. If None, max_mode_frequency is set based on the input timestep: max_mode_frequency = 1/(5*full_timestep)

• full_timestep (float, optional, by default .0001) – The stepper advances the dynamics by a time interval dt=full_timestep. Applies to both linear and non-linear stepper. A GS equilibrium is calculated every full_timestep. Note that this input is overridden by ‘automatic_timestep’ if the latter is not set to False.

• max_internal_timestep (float, optional, by default .0001) – Each time advancement of one full_timestep is divided in several sub-steps, with size of, at most, max_internal_timestep. Such sub_step is used to advance the circuit equations (under the assumption of constant applied voltage during the full_timestep). Note that this input is overridden by ‘automatic_timestep’ if the latter is not set to False.

• plasma_resistivity (float, optional, by default 1e-6) – Resistivity of the plasma. Plasma resistance values for each of the domain grid points are 2*np.pi*plasma_resistivity*eq.R/(dR*dZ) where dR*dZ is the area of the domain element.

• plasma_norm_factor (float, optional, by default 1000) – The plasma current is re-normalised by this factor, to bring to a value more akin to those of the metal currents

• blend_hatJ (float, optional, by default 0) – optional coefficient which enables use a blended version of the normalised plasma current distribution when contracting the plasma lumped circuit eq. from the left. The blend combines the current distribution at time t with (a guess for) the one at time t+dt.

• dIydI (np.array of size (np.sum(plasma_domain_mask), n_metal_modes+1), optional) – dIydI_(i,j) = d(Iy_i)/d(I_j) This is the jacobian of the plasma current distribution Iy with respect to all independent metal currents (both active and vessel modes) and to the total plasma current This is provided if known, otherwise calculated here at the linearization eq

• automatic_timestep ((float, float) or False, optional, by default False) – If not False, this overrides inputs full_timestep and max_internal_timestep: the timescales of the linearised problem are used to set the size of the timestep. The input eq and profiles are used to calculate the fastest growthrate, t_growthrate, henceforth, full_timestep = automatic_timestep[0]*t_growthrate max_internal_timestep = automatic_timestep[1]*full_timestep

• mode_removal (bool, optional, by default True) – It True, vessel normal modes are dropped after dIydI is calculated Modes that couple with the plasma less than min_dIy_dI than the strongest mode, are dropped. This criterion is applied based on the actual dIydI, calculated on GS solutions.

• linearize (bool, optional, by default True) – Whether to set up the linearization of the evolutive problem

• fix_n_vessel_modes (int) – If -1, modes are selected based on max_mode_frequency, threshold_dIy_dI and min_dIy_dI. If a non-negative integer, the number of vessel modes is fixed accordingly. max_mode_frequency, threshold_dIy_dI and min_dIy_dI are not used.

• threshold_dIy_dI (float) – Threshold value to drop vessel modes. Modes that couple with the plasma more than min_dIy_dI than the strongest mode, are included. This criterion is applied based on dIydI_noGS.

• min_dIy_dI (float) – Threshold value to drop vessel modes. Modes that couple with the plasma less than min_dIy_dI than the strongest mode, are dropped. This criterion is applied based on dIydI_noGS.

• custom_coil_resist (np.array) – 1d array of resistance values for all machine conducting elements, including both active coils and passive structures If None, the values calculated by default in tokamak will be sourced and used.

• custom_self_ind (np.array) – 2d matrix of mutual inductances between all pairs of machine conducting elements, including both active coils and passive structures If None, the values calculated by default in tokamak will be sourced and used.

assign_currents(currents_vec, eq, profiles)

Assigns current values as in input currents_vec to eq.tokamak and plasma. The input eq and profiles are modified accordingly. The format of the input currents aligns with self.currents_vec: (active coil currents, selected vessel normal modes currents, total plasma current/plasma_norm_factor)

Parameters:

• currents_vec (np.array) – Input current values to be assigned, in terms of mode currents

• eq (FreeGSNKE equilibrium Object) – Equilibrium object to be modified.

• profiles (FreeGSNKE profiles Object) – profiles object to be modified.

assign_currents_solve_GS(currents_vec, rtol_NK)

Assigns current values as in input currents_vec to private self.eq2 and self.profiles2. Static GS problem is accordingly solved, which finds the associated plasma flux and current distribution.

Parameters:

• currents_vec (np.array) – Input current values to be assigned. Format as in self.assign_currents.

• rtol_NK (float) – Relative tolerance to be used in the static GS problem.

build_current_vec(eq, profiles)

Builds the vector of currents in which the dynamics is actually solved, self.currents_vec This contains, in the order: (active coil currents, selected vessel normal modes currents, total plasma current/plasma_norm_factor)

Parameters:

• profiles (FreeGSNKE profiles Object) – profiles function of the initial equilibrium. Used to extract the value of the total plasma current.

• eq (FreeGSNKE equilibrium Object) – Initial equilibrium. Used to extract the value of all metal currents.

build_dIydI_j(j, rtol_NK)

Computes the term d(Iy)/dI_j of the Jacobian as a finite difference derivative, using the value of delta(I_j) inferred earlier by self.prepare_build_dIydI_j.

Parameters:

• j (int) – Index identifying the current to be varied. Indexes as in self.currents_vec.

• rtol_NK (float) – Relative tolerance to be used in the static GS problems.

Returns:

dIydIj – This is a 1d vector including all grid points in reduced domain, as from plasma_domain_mask.

Return type:

np.array finite difference derivative d(Iy)/dI_j.

build_dIydI_noGS(force_core_mask_linearization, starting_dI, core_mask, verbose)

Calculates a first estimate of the norm of dIy/dI for each mode without solving GS, i.e. using only the modified psi_tokamak. This is used in the mode selection for a first sifting of the modes. If force_core_mask_linearization is True, then alters self.starting_dI and self.approved_target_dIy so that the core mask is preserved.

Parameters:

• core_mask (np.array) – core mask of the reference equilibrium plasma

• force_core_mask_linearization (bool) – whether finite difference calculations should all be based on plasmas with the exact same core region

build_linearization(eq, profiles, dIydI, target_relative_tolerance_linearization, force_core_mask_linearization, verbose)

Builds the Jacobian d(Iy)/dI to set up the solver of the linearised problem.

Parameters:

• eq (FreeGS4E equilibrium Object) – Equilibrium around which to linearise.

• profiles (FreeGS4E profiles Object) – profiles properties of the equilibrium around which to linearise.

• dIydI (np.array) – input Jacobian, enter where available, otherwise this will be calculated here

• rtol_NK (float) – Relative tolerance to be used in the static GS problems.

• target_dIy (float, by default 0.001.) – Target relative value for the norm of delta(I_y), on which the finite difference derivative is calculated.

calc_lumped_plasma_resistance(norm_red_Iy0, norm_red_Iy1)

Uses the plasma resistance matrix R_yy to calculate the lumped plasma resistance, by contracting this with the vectors norm_red_Iy0, norm_red_Iy0. These should be normalised plasma current distribution vectors.

Parameters:

• norm_red_Iy0 (np.array) – Normalised plasma current distribution. This vector should sum to 1.

• norm_red_Iy1 (np.array) – Normalised plasma current distribution. This vector should sum to 1.

Returns:

Lumped resistance of the plasma.

Return type:

float

calculate_hatIy(trial_currents, plasma_psi)

Finds the normalised plasma current distribution corresponding to the combination of the input current values and plasma flux. Note that this does not assume that current values and plasma flux together are a solution of GS.

Parameters:

• trial_currents (np.array) – Vector of current values. Same format as self.currents_vec.

• plasma_psi (np.array) – Plasma flux values on full domain of shape (eq.nx, eq.ny), 2d.

Returns:

Normalised plasma current distribution. 1d vector on the reduced plasma domain.

Return type:

np.array

calculate_hatIy_GS(trial_currents, rtol_NK, record_for_updates=False)

Finds the normalised plasma current distribution corresponding to the combination of the input current values by solving the static GS problem.

Parameters:

• trial_currents (np.array) – Vector of current values. Same format as self.currents_vec.

• rtol_NK (float) – Relative tolerance to be used in the static GS problem.

Returns:

Normalised plasma current distribution. 1d vector on the reduced plasma domain.

Return type:

np.array

calculate_rel_tolerance_GS(trial_plasma_psi, a_res_GS=None)

Calculates how the residual in the plasma flux due to the static GS problem compares to the change in the plasma flux itself due to the dynamics, i.e. to the difference psi(t+dt) - psi(t) The relative residual is used to quantify the relative convergence of the stepper. It accesses self.trial_plasma_psi, self.eq1.plasma_psi, self.tokamak_psi

Parameters:

• trial_plasma_psi (ndarray) – psi(t+dt)

• a_res_GS (ndarray) – The residual of the static GS problem at t+dt

Returns:

Relative plasma flux residual.

Return type:

float

calculate_rel_tolerance_currents(current_residual, curr_eps)

Calculates how the current_residual in input compares to the step in the currents themselves, i.e. to the difference currents(t+dt) - currents(t-dt) This relative residual is used to quantify the relative convergence of the stepper. It accesses self.trial_currents and self.currents_vec_m1.

Parameters:

• current_residual (np.array) – Residual in current values. Same format as self.currents_vec.

• curr_eps (float) – Min value of the current step. Avoids divergence when dividing by the step in the currents.

Returns:

Relative current residual. Same format as self.currents_vec.

Return type:

np.array

check_and_change_plasma_resistivity(plasma_resistivity, relative_threshold_difference=0.01)

Checks if the plasma resistivity is different and resets it.

Parameters:

plasma_resistivity (float) – Resistivity of the plasma. Plasma resistance values for each of the domain grid points are 2*np.pi*plasma_resistivity*eq.R/(dR*dZ) where dR*dZ is the area of the domain element.

check_and_change_profiles(profiles_parameters=None)

Checks if new input parameters are provided for the profiles at t+dt. If so, it actions the necessary changes.

Parameters:

profiles_parameters (None or dictionary) – Set to None when the profiles parameter are left unchanged. Dictionary otherwise. See ‘get_profiles_values’ for dictionary structure.

currents_from_hatIy(hatIy1, active_voltage_vec)

Uses a guess for the normalised plasma current distribution at time t+dt to obtain all current values at time t+dt, through the ‘simplified’ circuit equations.

Parameters:

• hatIy1 (np.array) – Guess for the normalised plasma current distribution at time t+dt. Should be a vector that sums to 1. Reduced plasma domain only.

• active_voltage_vec (np.array) – Vector of active voltages for the active coils, applied between t and t+dt.

Returns:

Current values at time t+dt. Same format as self.currents_vec.

Return type:

np.array

get_profiles_values(profiles)

Extracts profiles properties.

Parameters:

profiles (FreeGS4E profiles Object) – profiles function of the initial equilibrium.

get_vessel_currents(eq)

Uses the input equilibrium to extract values for all metal currents, including active coils and vessel passive structures. These are stored in self.vessel_currents_vec

Parameters:

eq (FreeGSNKE equilibrium Object) – Initial equilibrium. eq.tokamak is used to extract current values.

hatIy1_iterative_cycle(hatIy1, active_voltage_vec, rtol_NK)

Uses a guess for the normalised plasma current distribution at time t+dt to obtain all current values at time t+dt through the circuit equations. Static GS is then solved for the same currents, which results in calculating the ‘iterated’ plasma flux and plasma current distribution at time t+dt (stored in the private self.eq2 and self.profiles2).

Parameters:

• hatIy1 (np.array) – Guess for the normalised plasma current distribution at time t+dt. Should be a vector that sums to 1. Reduced plasma domain only.

• active_voltage_vec (np.array) – Vector of active voltages for the active coils, applied between t and t+dt.

• rtol_NK (float) – Relative tolerance to be used in the static GS problem.

initialize_from_ICs(eq, profiles, target_relative_tolerance_linearization=1e-07, dIydI=None, force_core_mask_linearization=False, verbose=False)

Uses the input equilibrium and profiles as initial conditions and prepares to solve for their dynamics. If needed, sets the the linearised solver by calculating the Jacobian dIy/dI.

Parameters:

• eq (FreeGS4E equilibrium Object) – Initial equilibrium. This assigns all initial metal currents.

• profiles (FreeGS4E profiles Object) – profiles function of the initial equilibrium. This assigns the initial total plasma current.

• rtol_NK (float) – Relative tolerance to be used in the static GS problems in the initialization and when calculating the Jacobian dIy/dI to set up the linearised problem. This does not set the tolerance of the static GS solves used by the dynamic solver, which is set through the stepper itself.

• dIydI (np.array of size (np.sum(plasma_domain_mask), n_metal_modes+1), optional) – dIydI_(i,j) = d(Iy_i)/d(I_j) This is the jacobian of the plasma current distribution with respect to all independent metal currents (both active and vessel modes) and to the total plasma current If not provided, this is calculated based on the properties of the provided equilibrium.

• verbose (bool) – Whether to allow progress printouts

make_blended_hatIy_(hatIy1, blend)

Averages the normalised plasma current distributions at time t and (a guess for the one at) at time t+dt to better contract the system of plasma circuit eqs.

Parameters:

hatIy1 (np.array) – Guess for the normalised plasma current distributions at time t+dt. Should be a vector that sums to 1. Reduced plasma domain only.

nlstepper(active_voltage_vec, profiles_parameters=None, plasma_resistivity=None, target_relative_tol_currents=0.005, target_relative_tol_GS=0.005, working_relative_tol_GS=0.001, target_relative_unexplained_residual=0.5, max_n_directions=3, step_size_psi=2.0, step_size_curr=0.8, scaling_with_n=0, blend_GS=0.5, curr_eps=1e-05, max_no_NK_psi=5.0, clip=5, verbose=0, linear_only=False)

This is the main stepper function. If linear_only = True, this advances the linearised dynamic problem. If linear_only = False, a solution of the full non-linear problem is seeked using a combination of NK methods. When a solution has been found, time is advanced by self.dt_step, the new values of all extensive currents are recorded in self.currents_vec and new equilibrium and profiles properties in self.eq1 and self.profiles1.

The solver’s algorithm proceeds like below: 1) solve linearised problem to obtain an initial guess of the currents and solve the associated static GS problem, assign such trial_plasma_psi and trial_currents (including the resulting tokamak_psi); 2) if pair [trial_plasma_psi, tokamak_psi] fails static GS tolerance check, update trial_plasma_psi bringing it closer to the actual GS solution; 3) at fixed trial_currents (and associated tokamak_psi) update trial_plasma_psi using NK solver for the associated root problem; 4) at fixed trial_plasma_psi, update trial_currents (and associated tokamak_psi) using NK solver for the associated root problem; 5) if convergence on the current residuals is not reached or static GS tolerance check fails, restart from point 2; 6) the pair [trial_currents, trial_plasma_psi] solves the nonlinear dynamic problem, assign values to self.currents_vec, self.eq1 and self.profiles1.

Parameters:

• active_voltage_vec (np.array) – Vector of active voltages for the active coils, applied between t and t+dt.

• profiles_parameters (Dictionary) – Set to None when the profiles parameters are left unchanged. Otherwise, dictionary containing the relevant profiles parameters for the profiles object on which the evolution is calculated See ‘get_profiles_values’ for dictionary structure.

• target_relative_tol_currents (float, optional, by default .005) – Relative tolerance in the currents required for convergence of the dynamic problem. This is calculated with respect to the change in the currents themselves due to the dynamical evolution: residual/(currents(t+dt) - currents(t-dt))

• target_relative_tol_GS (float, optional, by default .005) – Relative tolerance in the plasma flux for the static GS problem required for convergence. This is calculated with respect to the change in the flux itself due to the dynamical evolution: residual/delta(psi(t+dt) - psi(t))

• working_relative_tol_GS (float, optional, by default .001) – Tolerance used when solving all static GS problems while executing the step, also expressed in relative terms as target_relative_tol_GS. Note this value needs to be smaller than target_relative_tol_GS to allow for convergence.

• target_relative_unexplained_residual (float, optional, by default .5) – Used in the NK solvers. Inclusion of additional Krylov basis vectors is stopped if the fraction of the residual (linearly) canceled is > 1-target_relative_unexplained_residual.

• max_n_directions (int, optional, by default 3) – Used in the NK solvers. Inclusion of additional Krylov basis vectors is stopped if max_n_directions have already been included.

• step_size_psi (float, optional, by default 2.) – Used by the NK solver applied to the root problem in the plasma flux. l2 norm of proposed step for the finite difference calculation, in units of the residual.

• step_size_curr (float, optional, by default .8) – Used by the NK solver applied to the root problem in the currents. l2 norm of proposed step for the finite difference calculation, in units of the residual.

• scaling_with_n (int, optional, by default 0) – Used in the NK solvers. Allows to further scale dx candidate steps by factor (1 + self.n_it)**scaling_with_n

• max_no_NK_psi (float, optional, by default 5.) – Execution of NK update on psi for the dynamic problem is triggered when relative_psi_residual > max_no_NK_psi * target_relative_tol_GS where the psi residual is calculated as for target_relative_tol_GS

• blend_GS (float, optional, by default .5) – Blend coefficient used in trial_plasma_psi updates at step 2 of the algorithm above. Should be between 0 and 1.

• curr_eps (float, optional, by default 1e-5) – Regulariser used in calculating the relative convergence on the currents. Avoids divergence when dividing by the advancement in the currents. Min value of the current per step.

• clip (float, optional, by default 5) – Used in the NK solvers. Maximum step size for each accepted basis vector, in units of the exploratory step.

• verbose (int, optional, by default 0) – Printouts of convergence process. Use 1 for printouts with details on each NK cycle. Use 2 for printouts with deeper intermediate details.

• linear_only (bool, optional, by default False) – If linear_only = True the solution of the linearised problem is accepted. If linear_only = False, the convergence criteria are used and a solution of the full nonlinear problem is seeked.

prepare_build_dIydI_j(j, rtol_NK, target_dIy, starting_dI, GS=True)

Prepares to compute the term d(Iy)/dI_j of the Jacobian by inferring the value of delta(I_j) corresponding to a change delta(I_y) with norm(delta(I_y))=target_dIy.

Parameters:

• j (int) – Index identifying the current to be varied. Indexes as in self.currents_vec.

• rtol_NK (float) – Relative tolerance to be used in the static GS problems.

• target_dIy (float) – Target value for the norm of delta(I_y), on which th finite difference derivative is calculated.

• starting_dI (float) – Initial value to be used as delta(I_j) to infer the slope of norm(delta(I_y))/delta(I_j).

• min_curr (float, optional, by default 1e-4) – If inferred current value is below min_curr, clip to min_curr.

• max_curr (int, optional, by default 300) – If inferred current value is above min_curr, clip to max_curr.

remove_modes(eq, selected_modes_mask)

It actions the removal of the unselected normal modes. Given a setup with a set of normal modes and a resulting size of the vector self.currents_vec, modes that are not selected in the input mask are removed and the circuit equations updated accordingly. The dimensionality of the vector self.currents_vec is reduced.

Parameters:

• eq (class) – FreeGSNKE equilibrium object.

• selected_modes_mask (np.array of bool values,) – shape(selected_modes_mask) = shape(self.currents_vec) at the time of calling the function indexes corresponding to True are kept, indexes corresponding to False are dropped

reset_plasma_resistivity(plasma_resistivity)

Function to reset the resistivity of the plasma.

Parameters:

plasma_resistivity (float) – Resistivity of the plasma. Plasma resistance values for each of the domain grid points are 2*np.pi*plasma_resistivity*eq.R/(dR*dZ) where dR*dZ is the area of the domain element.

reset_timestep(full_timestep, max_internal_timestep)

Allows for a reset of the timesteps.

Parameters:

• full_timestep (float) – The stepper advances the dynamics by a time interval dt=full_timestep. Applies to both linear and non-linear stepper. A GS equilibrium is calculated every full_timestep.

• max_internal_timestep (float) – Each time advancement of one full_timestep is divided in several sub-steps, with size of, at most, max_internal_timestep. Such sub_step is used to advance the circuit equations

set_linear_solution(active_voltage_vec, d_profiles_pars_dt=None)

Uses the solver of the linearised problem to set up an initial guess for the nonlinear solver for the currents at time t+dt. Uses self.currents_vec as I(t). Solves GS at time t+dt using the currents derived from the linearised dynamics.

Parameters:

• active_voltage_vec (np.array) – Vector of external voltage applied to the active coils during the timestep.

• d_profiles_pars_dt (None) – Does not currently use d_profiles_pars_dt The evolution of the profiles coefficient is not accounted by the linearised dynamical evolution.

set_plasma_resistivity(plasma_resistivity)

Function to set the resistivity of the plasma. self.plasma_resistance_1d is the diagonal of the matrix R_yy, the plasma resistance matrix. Note that it only spans the grid points in the reduced domain, as from plasma_domain_mask.

Parameters:

plasma_resistivity (float) – Resistivity of the plasma. Plasma resistance values for each of the domain grid points are 2*np.pi*plasma_resistivity*eq.R/(dR*dZ) where dR*dZ is the area of the domain element.

step_complete_assign(working_relative_tol_GS, from_linear=False)

This function completes the advancement by dt. The time-evolved currents as obtained by the stepper (self.trial_currents) are recorded in self.currents_vec and assigned to the equilibrium self.eq1. The time-evolved equilibrium properties (i.e. the plasma flux self.trial_plasma_psi and resulting current distribution) are recorded in self.eq1 and self.profiles1.

Parameters:

• working_relative_tol_GS (float) – The relative tolerance of the GS solver used to solve the dynamics is set to a fraction of the change in the plasma flux associated to the timestep itself. The fraction is set through working_relative_tol_GS.

• from_linear (bool, optional, by default False) – If the stepper is only solving the linearised problem, use from_linear=True. This acellerates calculations by reducing the number of static GS solve calls.

unstable_mode_deformations(starting_dI=50, rtol_NK=1e-07, target_dIy=0.002)

Applies the unstable mode m to calculate (dR/dIm, dZ/dIm) where R and Z are the current-averaged coords of the plasma (i.e. eq.Rcurrent and eq.Zcurrent) and Im is the magnitude of the current in the unstable mode. A map of the current distribution after Im is applied is recorded, together with the map having the same R and Z obtained as a rigid displament of the original eq.

Parameters:

• rtol_NK (float) – Relative tolerance to be used in the static GS problems.

• target_dIy (float) – Target value for the norm of delta(I_y), on which the finite difference derivative is calculated.

• starting_dI (float) – Initial value to be used as delta(I_j) to infer the slope of norm(delta(I_y))/delta(Im).