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)