freegsnke.GSstaticsolver module

Applies the Newton Krylov solver Object to the static GS problem. Implements both forward and inverse GS solvers.

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.GSstaticsolver.NKGSsolver(eq)[source]

Bases: object

Solver for the non-linear forward Grad Shafranov (GS) static problem. Here, the GS problem is written as a root problem in the plasma flux psi. This root problem is passed to and solved by the NewtonKrylov solver itself, class nk_solver.

The solution domain is set at instantiation time, through the input FreeGSNKE equilibrium object.

The non-linear solvers are called using the ‘forward_solve’, ‘inverse_solve’ or generic ‘solve’ methods.

F_function(plasma_psi, tokamak_psi, profiles)[source]

Residual of the nonlinear Grad Shafranov equation written as a root problem F(plasma_psi) equiv [delta* - J](plasma_psi) The plasma_psi that solves the Grad Shafranov problem satisfies F(plasma_psi) = [delta* - J](plasma_psi) = 0

Parameters:

plasma_psi (np.array of size eq.nx*eq.ny) – magnetic flux due to the plasma
tokamak_psi (np.array of size eq.nx*eq.ny) – magnetic flux due to the tokamak alone, including all metal currents, in both active coils and passive structures
profiles (freegs4e profile object) – profile object describing target plasma properties, used to calculate current density jtor

Returns:

residual – residual of the GS equation

Return type:

np.array of size eq.nx*eq.ny

__init__(eq)[source]

Instantiates the solver object. Based on the domain grid of the input equilibrium object, it prepares - the linear solver ‘self.linear_GS_solver’ - the response matrix of boundary grid points ‘self.greens_boundary’

Parameters:: eq (a FreeGSNKE equilibrium object.) – The domain grid defined by (eq.R, eq.Z) is the solution domain adopted for the GS problems. Calls to the nonlinear solver will use the grid domain set at instantiation time. Re-instantiation is necessary in order to change the propertes of either grid or domain.

forward_solve(eq, profiles, target_relative_tolerance, max_solving_iterations=100, Picard_handover=0.11, step_size=2.5, scaling_with_n=-1.0, target_relative_unexplained_residual=0.2, max_n_directions=16, max_rel_update_size=0.2, clip=10, force_up_down_symmetric=False, verbose=False, suppress=False)[source]

The method that actually solves the forward static GS problem.

A forward problem is specified by the 2 FreeGSNKE objects eq and profiles. The first specifies the metal currents (throught eq.tokamak) and the second specifies the desired plasma properties (i.e. plasma current and profile functions).

The plasma_psi which solves the given GS problem is assigned to the input eq, and can be found at eq.plasma_psi.

Parameters:

eq (FreeGSNKE equilibrium object) – Used to extract the assigned metal currents, which in turn are used to calculate the according self.tokamak_psi
profiles (FreeGSNKE profile object) – Specifies the target properties of the plasma. These are used to calculate Jtor(psi)
target_relative_tolerance (float) – NK iterations are interrupted when this criterion is satisfied. Relative convergence for the residual F(plasma_psi)
max_solving_iterations (int) – NK iterations are interrupted when this limit is surpassed
Picard_handover (float) – Value of relative tolerance above which a Picard iteration is performed instead of a NK call
step_size (float) – l2 norm of proposed step, in units of the size of the residual R0
scaling_with_n (float) – allows to further scale the proposed steps as a function of the number of previous steps already attempted (1 + self.n_it)**scaling_with_n
target_relative_unexplained_residual (float between 0 and 1) – terminates internal iterations when the considered directions can (linearly) explain such a fraction of the initial residual R0
max_n_directions (int) – terminates iteration even though condition on explained residual is not met
max_rel_update_size (float) – maximum relative update, in norm, to plasma_psi. If larger than this, the norm of the update is reduced
clip (float) – maximum size of the update due to each explored direction, in units of exploratory step used to calculate the finite difference derivative
force_up_down_symmetric (bool) – If True, the solver is forced on up-down symmetric trial solutions
verbose (bool) – flag to allow progress printouts
suppress (bool) – flag to allow suppress all printouts

freeboundary(plasma_psi, tokamak_psi, profiles)[source]

Imposes boundary conditions on set of boundary points.

Parameters:

plasma_psi (np.array of size eq.nx*eq.ny) – magnetic flux due to the plasma
tokamak_psi (np.array of size eq.nx*eq.ny) – magnetic flux due to the tokamak alone, including all metal currents, in both active coils and passive structures
profiles (FreeGSNKE profile object) – profile object describing target plasma properties. Used to calculate current density jtor

get_rel_delta_psi(new_psi, previous_psi, profiles)[source]

Calculates the relative change between new_psi and previous_psi in the core region

Parameters:

new_psi (np.array) – Flattened flux function
previous_psi (np.array) – _descrFlattened flux functioniption_
profiles (freegsnke profile object) – Used to source the core mask

get_rel_delta_psit(delta_current, profiles, vgreen)[source]

Calculates the relative change to the tokamak_psi in the core region caused by the requested current change ‘delta_current’.

Parameters:

delta_current (np.array) – Vector of requested current changes.
profiles (freegsnke profile object) – Used to source the core mask
vgreen (np.array) – The green functions of the relevant coils. For example eq._vgreen

inverse_solve(eq, profiles, constrain, target_relative_tolerance, target_relative_psit_update=0.001, max_solving_iterations=100, max_iter_per_update=5, Picard_handover=0.15, step_size=2.5, scaling_with_n=-1.0, target_relative_unexplained_residual=0.3, max_n_directions=16, clip=10, max_rel_update_size=0.15, threshold_val=0.18, l2_reg=1e-09, forward_tolerance_increase=100, max_rel_psit=0.02, damping_factor=0.995, use_full_Jacobian=True, full_jacobian_handover=[1e-05, 0.007], l2_reg_fj=1e-08, force_up_down_symmetric=False, verbose=False, suppress=False)[source]

Inverse solver.

Both coil currents and plasma flux are sought. Needs a set of desired magnetic constraints.

An inverse problem is specified by the 2 FreeGSNKE objects, eq and profiles, plus a freegsnke “constrain” (or Inverse_optimizer) object. The first specifies the metal currents (throught eq.tokamak) The second specifies the desired plasma properties (i.e. plasma current and profile functions). The constrain object collects the desired magnetic constraints.

The plasma_psi which solves the given GS problem is assigned to the input eq, and can be found at eq.plasma_psi. The coil currents with satisfy the magnetic constraints are assigned to eq.tokamak

Parameters:

eq (FreeGSNKE equilibrium object) – Used to extract the assigned metal currents, which in turn are used to calculate the according self.tokamak_psi
profiles (FreeGSNKE profile object) – Specifies the target properties of the plasma. These are used to calculate Jtor(psi)
constrain (freegsnke inverse_optimizer object) – Specifies the coils available for control and the desired magnetic constraints
target_relative_tolerance (float) – The desired tolerance for the plasma flux. At fixed coil currents, this is the tolerance imposed to the GS problem. This has to be satisfied for the inverse problem to be considered solved.
target_relative_psit_update (float) – The relative update to psi_tokamak caused by the update in the control currents has to be lower than this target value for the inverse problem to be considered successfully solved.
max_solving_iterations (int) – Maximum number of solving iterations. The solver is interuupted when this is met.
max_iter_per_update (int) – Maximum number of interations allowed to the forward solver in each cycle of the inverse solver.
Picard_handover (float) – Value of relative tolerance above which a Picard iteration is performed instead of a NK call in the forward solve
step_size (float) – l2 norm of proposed step, in units of the size of the residual R0
scaling_with_n (float) – allows to further scale the proposed steps as a function of the number of previous steps already attempted (1 + self.n_it)**scaling_with_n
target_relative_unexplained_residual (float between 0 and 1) – terminates internal iterations when the considered directions can (linearly) explain such a fraction of the initial residual R0
max_n_directions (int) – terminates iteration even though condition on explained residual is not met
max_rel_update_size (float) – maximum relative update, in norm, to plasma_psi. If larger than this, the norm of the update is reduced
clip (float) – maximum size of the update due to each explored direction, in units of exploratory step used to calculate the finite difference derivative
l2_reg (either float or 1d np.array with len=self.n_control_coils) – The regularization factor applied when green functions are used as Jacobians
forward_tolerance_increase (float) – after coil currents are updated, the interleaved forward problems are requested to converge to a tolerance that is tighter by a factor forward_tolerance_increase with respect to the change in flux caused by the current updates over the plasma core
forward_tolerance_increase_factor (float) – iterations that do not result in improvement trigger an increase in the tolerance applied to the forward problems. The increase corresponds to a factor forward_tolerance_increase_factor
max_rel_psit (float) – The maximum relative change that the requested updates in the control currents are allowed to cause.
damping_factor (float) – This applies a damping of damping_factor**number_of_iterations to max_rel_psit, to encourage convergence
full_jacobian_handover (float) – When the forward problems achieve this tolerance level, self.optimize_currents is called instead of constrain.optimize_currents. This means that the actual Jacobians are used rather than the green functions.
l2_reg_fj (float) – The regularization factor applied when the full Jacobians are used.
verbose (bool) – flag to allow progress printouts
suppress (bool) – flag to allow suppress all printouts

optimize_currents(eq, profiles, constrain, target_relative_tolerance, relative_psit_size=0.001, l2_reg=1e-12, verbose=False)[source]

Calculates requested current changes for the coils available to control using the actual Jacobian rather than assuming the Jacobian is given by the green functions.

Parameters:

eq (freegsnke equilibrium object) – Used to source coil currents and plasma_psi
profiles (freegsnke equilibrium object) – Provides the plasma properties
constrain (freegnske inverse_optimizer object) – Provides information on the coils available for control
target_relative_tolerance (float) – Target tolerance applied to GS when building the Jacobian
relative_psit_size (float, optional) – Used to size the finite difference steps that define the Jacobian, by default 1e-3
l2_reg (float, optional) – The Tichonov regularization factor applied to the least sq problem, by default 1e-12
verbose (bool) – Print output

port_critical(eq, profiles)[source]

Transfers critical points and other useful info from profile to equilibrium object, after GS solution.

Parameters:

eq (FreeGSNKE equilibrium object) – Equilibrium on which to record values
profiles (FreeGSNKE profile object) – Profiles object which has been used to calculate Jtor.

relative_del_residual(res, psi)[source]

Calculates a normalised relative residual, based on the norm max(.) - min(.)

Parameters:

res (ndarray) – Residual
psi (ndarray) – plasma_psi

Returns:

Relative normalised residual, norm(plasma_psi)

Return type:

float, float

relative_norm_residual(res, psi)[source]

Calculates a normalised relative residual, based on linalg.norm

Parameters:

res (ndarray) – Residual
psi (ndarray) – plasma_psi

Returns:

Relative normalised residual

Return type:

float

solve(eq, profiles, constrain=None, target_relative_tolerance=1e-05, target_relative_psit_update=0.001, max_solving_iterations=100, max_iter_per_update=5, Picard_handover=0.1, step_size=2.5, scaling_with_n=-1.0, target_relative_unexplained_residual=0.3, max_n_directions=16, clip=10, max_rel_update_size=0.15, l2_reg=1e-09, forward_tolerance_increase=100, max_rel_psit=0.01, damping_factor=0.98, use_full_Jacobian=True, full_jacobian_handover=[1e-05, 0.007], l2_reg_fj=1e-08, force_up_down_symmetric=False, verbose=False, suppress=False)[source]

The method to solve the GS problems, both forward and inverse.

an inverse solve is specified by the ‘constrain’ input,

which includes the desired constraints on the configuration of magnetic flux. Please see inverse_solve for details. - a forward solve has constrain=None. Please see forward_solve for details.

Parameters:

eq (FreeGSNKE equilibrium object) – Used to extract the assigned metal currents, which in turn are used to calculate the according self.tokamak_psi
profiles (FreeGSNKE profile object) – Specifies the target properties of the plasma. These are used to calculate Jtor(psi)
constrain (freegsnke inverse_optimizer object) – Specifies the coils available for control and the desired magnetic constraints
target_relative_tolerance (float) – The desired tolerance for the plasma flux. At fixed coil currents, this is the tolerance imposed to the GS problem. This has to be satisfied for the inverse problem to be considered solved.
target_relative_psit_update (float) – The relative update to psi_tokamak caused by the update in the control currents has to be lower than this target value for the inverse problem to be considered successfully solved.
max_solving_iterations (int) – Maximum number of solving iterations. The solver is interuupted when this is met.
max_iter_per_update (int) – Maximum number of interations allowed to the forward solver in each cycle of the inverse solver.
Picard_handover (float) – Value of relative tolerance above which a Picard iteration is performed instead of a NK call in the forward solve
step_size (float) – l2 norm of proposed step, in units of the size of the residual R0
scaling_with_n (float) – allows to further scale the proposed steps as a function of the number of previous steps already attempted (1 + self.n_it)**scaling_with_n
target_relative_unexplained_residual (float between 0 and 1) – terminates internal iterations when the considered directions can (linearly) explain such a fraction of the initial residual R0
max_n_directions (int) – terminates iteration even though condition on explained residual is not met
max_rel_update_size (float) – maximum relative update, in norm, to plasma_psi. If larger than this, the norm of the update is reduced
clip (float) – maximum size of the update due to each explored direction, in units of exploratory step used to calculate the finite difference derivative
l2_reg (either float or 1d np.array with len=self.n_control_coils) – The regularization factor applied when green functions are used as Jacobians
forward_tolerance_increase (float) – after coil currents are updated, the interleaved forward problems are requested to converge to a tolerance that is tighter by a factor forward_tolerance_increase with respect to the change in flux caused by the current updates over the plasma core
forward_tolerance_increase_factor (float) – iterations that do not result in improvement trigger an increase in the tolerance applied to the forward problems. The increase corresponds to a factor forward_tolerance_increase_factor
max_rel_psit (float) – The maximum relative change that the requested updates in the control currents are allowed to cause.
damping_factor (float) – This applies a damping of damping_factor**number_of_iterations to max_rel_psit, to encourage convergence
full_jacobian_handover (float) – When the forward problems achieve this tolerance level, self.optimize_currents is called instead of constrain.optimize_currents. This means that the actual Jacobians are used rather than the green functions.
l2_reg_fj (float) – The regularization factor applied when the full Jacobians are used.
verbose (bool) – flag to allow progress printouts
suppress (bool) – flag to allow suppress all printouts