freegsnke.linear_solve module

Implements the object that advances the linearised dynamics.

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.linear_solve.linear_solver(eq, Lambdam1, Pm1, Rm1, Mey, plasma_norm_factor, plasma_resistance_1d, max_internal_timestep=0.0001, full_timestep=0.0001)[source]

Bases: object

Interface between the linearised system of circuit equations and the implicit-Euler time stepper. Calculates the linear growth rate and solves the linearised dynamical problem. It needs the Jacobian of the plasma current distribution with respect to all of the independent currents, dIy/dI.

__init__(eq, Lambdam1, Pm1, Rm1, Mey, plasma_norm_factor, plasma_resistance_1d, max_internal_timestep=0.0001, full_timestep=0.0001)[source]

Instantiates the linear_solver object, with inputs computed mostly within the circuit_equation_metals object. Based on the input plasma properties and coupling matrices, it prepares: - an instance of the implicit Euler solver implicit_euler_solver() - internal time-stepper for the implicit-Euler

Parameters:
  • eq (class) – FreeGSNKE equilibrium object.

  • Lambdam1 (np.array) – State matrix of the circuit equations for the metal in normal mode form: P is the identity on the active coils and diagonalises the isolated dynamics of the passive coils, R^{-1/2}L_{passive}R^{-1/2}

  • Pm1 (np.array) – change of basis matrix, as defined above, to the power of -1

  • Rm1 (np.array) – matrix of all metal resitances to the power of -1. Diagonal.

  • Mey (np.array) – matrix of inductance values between grid points in the reduced plasma domain and all metal coils (active coils and passive-structure filaments) Calculated by the metal_currents object

  • plasma_norm_factor (float) – an overall factor to work with a rescaled plasma current, so that it’s within a comparable range with metal currents

  • max_internal_timestep (float) – internal integration timestep of the implicit-Euler solver, to be used as substeps over the <<full_timestep>> interval

  • full_timestep (float) – full timestep requested to the implicit-Euler solver

build_Mmatrix()[source]

Initialises the pseudo-inductance matrix of the problem Mdot(x) + Rx = forcing using the linearisation Jacobian.

Lambda^-1 + P^-1R^-1Mey A P^-1R^-1Mey B

M = M0 + dM = ( )

J(Myy A + MyeP)/Rp J Myy B/Rp

This also builds the forcing:

P^-1R^-1 Voltage P^-1R^-1Mey

forcing = ( ) - ( ) C dot{theta}

0 J Myy/Rp

where A = dIy/dId

B = dIy/dIp C = dIy/plasmapars

Here we take C=0, that is the linearised dynamics does not account for evolving plasma parameters atm.

Parameters:
  • explicitly (None given)

  • attributes. (they are all given by the object)

calculate_linear_growth_rate()[source]

Looks into the eigenvecotrs of the “M” matrix to find the negative singular values, which correspond to the growth rates of instabilities.

Parameters:

attributes (parameters are passed in as object)

calculate_pseudo_rigid_projections(dRZdI)[source]

Projects the unstable modes on the vectors of currents which best isolate an R or a Z movement of the plasma

Parameters:

dRZdI (np.array) – Jacobian of Rcurrent and Zcurrent shifts wrt the modes, as calculated in nonlinear_solve

Returns:

proj[i,0] is the scalar product of the unstable mode i on the vector of modes resulting in an Rcurrent shift proj[i,1] is the scalar product of the unstable mode i on the vector of modes resulting in an Zcurrent shift

Return type:

np.array

calculate_stability_margin()[source]

Here we calculate the stability margin parameter from:

https://iopscience.iop.org/article/10.1088/0029-5515/45/8/021

Parameters:

attributes (parameters are passed in as object)

reset_plasma_resistivity(plasma_resistance_1d)[source]

Resets the value of the plasma resistivity, throught the vector of ‘geometric restistances’ in the plasma domain

Parameters:

plasma_resistance_1d (ndarray) – Vector of 2pi resistivity R/dA for all domain grid points in the reduced plasma domain

reset_timesteps(max_internal_timestep, full_timestep)[source]

Resets the integration timesteps, calling self.solver.set_timesteps

Parameters:
  • max_internal_timestep (float) – integration substep of the ciruit equation, calling an implicit-Euler solver

  • full_timestep (float) – integration timestep of the circuit equation

set_linearization_point(dIydI, hatIy0, Myy_hatIy0)[source]

Initialises an implicit-Euler solver with the appropriate matrices for the linearised dynamic problem.

Parameters:
  • dIydI (np.array) – partial derivatives of plasma-cell currents on the reduced plasma domain with respect to all intependent metal currents (active coil currents, vessel normal modes, total plasma current divided by plasma_norm_factor). These would typically come from having solved the forward Grad-Shafranov problem. Finite difference Jacobian. Calculated by the nl_solver object

  • hatIy0 (np.array) – Plasma current distribution on the reduced plasma domain (1d) of the equilibrium around which the dynamics is linearised. This is normalised by the total plasma current, so that this vector sums to 1.

  • Myy_hatIy0 (np.array) – The matrix product np.dot(Myy, hatIy0) in the same reduced domain as hatIy0 This is provided by Myy_handler

stepper(It, active_voltage_vec, d_profiles_pars_dt=None)[source]

Executes the time advancement. Uses the implicit_euler instance.

Parameters:
  • It (np.array) – vector of all independent currents that are solved for by the linearides problem, in terms of normal modes: (active currents, vessel normal modes, total plasma current divided by normalisation factor)

  • active_voltage_vec (np.array) – voltages applied to the active coils

  • d_profiles_pars_dt (np.array) – time derivative of the profile parameters, not used atm

  • attributes (other parameters are passed in as object)