Coverage for python/gsfit/database_readers/st40_mdsplus_with_digital_filter/setup_plasma.py: 0%

69 statements  

« prev     ^ index     » next       coverage.py v7.16.2, created at 2026-10-14 15:25 +0000

1import typing 

2from typing import TYPE_CHECKING 

3 

4import gsfit_rs 

5import numpy as np 

6import numpy.typing as npt 

7from gsfit_rs import Plasma 

8from st40_database import GetData 

9 

10from ...tensioned_cubic_splines_regularisations import make_tensioned_cubic_b_spline_regularisations 

11 

12if TYPE_CHECKING: 

13 from . import DatabaseReader 

14 

15 

16def setup_plasma( 

17 self: "DatabaseReader", 

18 pulseNo: int, 

19 settings: dict[str, typing.Any], 

20 times_to_reconstruct: npt.NDArray[np.float64], 

21) -> Plasma: 

22 """ 

23 This method initialises the Rust `Plasma` class. 

24 

25 :param pulseNo: Pulse number, used to read from the database 

26 :param settings: Dictionary containing the JSON settings read from the `settings` directory 

27 :param times_to_reconstruct: Times the equilibrium will be solved at [second] 

28 

29 **This method is specific to ST40's experimental MDSplus database.** 

30 

31 See `python/gsfit/database_readers/interface.py` for more details on how a new database_reader should be implemented. 

32 """ 

33 

34 # Initial plasma conditions 

35 initial_guess_ip = settings["GSFIT_code_settings.json"]["initial_guess"]["ip"] 

36 initial_guess_cur_r = settings["GSFIT_code_settings.json"]["initial_guess"]["cur_r"] 

37 initial_guess_cur_z = settings["GSFIT_code_settings.json"]["initial_guess"]["cur_z"] 

38 initial_guess_minor_radius = settings["GSFIT_code_settings.json"]["initial_guess"]["minor_radius"] 

39 initial_guess_elongation = settings["GSFIT_code_settings.json"]["initial_guess"]["elongation"] 

40 

41 # Numerical settings the Grad-Shafranov solve is run with 

42 n_iter_max = settings["GSFIT_code_settings.json"]["numerics"]["n_iter_max"] 

43 n_iter_min = settings["GSFIT_code_settings.json"]["numerics"]["n_iter_min"] 

44 n_iter_no_vertical_feedback = settings["GSFIT_code_settings.json"]["numerics"]["n_iter_no_vertical_feedback"] 

45 grad_shafranov_deviation_tolerance = settings["GSFIT_code_settings.json"]["numerics"]["grad_shafranov_deviation_tolerance"] 

46 use_anderson_mixing = settings["GSFIT_code_settings.json"]["numerics"]["anderson_mixing"]["use"] 

47 anderson_mixing_from_previous_iter = settings["GSFIT_code_settings.json"]["numerics"]["anderson_mixing"]["mixing_from_previous_iter"] 

48 # The pressure at the plasma boundary. Integrating `p'` leaves the pressure defined up to a constant, which 

49 # GSFit historically fixed by demanding zero pressure at the boundary. Setting `free` makes it a degree of 

50 # freedom of the fit instead. `.get` is used so that older settings directories keep working 

51 pressure_edge_settings = settings["source_function_p_prime.json"].get( 

52 "pressure_edge", settings["GSFIT_code_settings.json"]["numerics"].get("pressure_edge", {}) 

53 ) 

54 pressure_edge_free = pressure_edge_settings.get("free", False) 

55 pressure_edge_regularisation_weight = pressure_edge_settings.get("regularisation_weight", 1.0e-4) 

56 

57 # Set the source functions types 

58 p_prime_source_function: gsfit_rs.EfitPolynomial | gsfit_rs.TensionedCubicBSpline 

59 ff_prime_source_function: gsfit_rs.EfitPolynomial | gsfit_rs.TensionedCubicBSpline 

60 

61 if settings["source_function_p_prime.json"]["method"] == "efit_polynomial": 

62 n_dof = settings["source_function_p_prime.json"]["efit_polynomial"]["n_dof"] 

63 regularisations = np.array(settings["source_function_p_prime.json"]["efit_polynomial"]["regularizations"]) 

64 # If `regularisations` is [[]] in the json file, will be interpreted by numpy as having size (1, 0). 

65 # Which would be interpreted as (n_regularisations, n_dof). So it would cause an error 

66 if regularisations.shape == (1, 0): 

67 regularisations = np.zeros((0, n_dof), dtype=np.float64) 

68 p_prime_source_function = gsfit_rs.EfitPolynomial(n_dof, regularisations) 

69 elif settings["source_function_p_prime.json"]["method"] == "tensioned_cubic_b_spline": 

70 p_prime_tensioned_cubic_b_spline_settings = settings["source_function_p_prime.json"]["tensioned_cubic_b_spline"] 

71 interior_knots = np.array(p_prime_tensioned_cubic_b_spline_settings["interior_knots"]) 

72 interval_tensions = np.array(p_prime_tensioned_cubic_b_spline_settings["interval_tensions"]) 

73 n_dof = len(interior_knots) + 4 

74 if p_prime_tensioned_cubic_b_spline_settings.get("make_cubic_spline_regularisations", False): 

75 # Build the regularisations from settings (see `tensioned_cubic_splines_regularisations.py`) 

76 regularisations = make_tensioned_cubic_b_spline_regularisations( 

77 interior_knots, 

78 interval_tensions, 

79 p_prime_tensioned_cubic_b_spline_settings["regularisation_builder"], 

80 ) 

81 else: 

82 regularisations = np.array(p_prime_tensioned_cubic_b_spline_settings["regularizations"]) 

83 # If `regularisations` is [[]] in the json file, will be interpreted by numpy as having size (1, 0). 

84 # Which would be interpreted as (n_regularisations, n_dof). So it would cause an error 

85 if regularisations.shape == (1, 0): 

86 regularisations = np.zeros((0, n_dof), dtype=np.float64) 

87 p_prime_source_function = gsfit_rs.TensionedCubicBSpline(regularisations, interior_knots, interval_tensions) 

88 else: 

89 raise ValueError(f"Unknown method for p_prime source function: {settings['source_function_p_prime.json']['method']}") 

90 

91 if settings["source_function_ff_prime.json"]["method"] == "efit_polynomial": 

92 n_dof = settings["source_function_ff_prime.json"]["efit_polynomial"]["n_dof"] 

93 regularisations = np.array(settings["source_function_ff_prime.json"]["efit_polynomial"]["regularizations"]) 

94 # If `regularisations` is [[]] in the json file, will be interpreted by numpy as having size (1, 0). 

95 # Which would be interpreted as (n_regularisations, n_dof). So it would cause an error 

96 if regularisations.shape == (1, 0): 

97 regularisations = np.zeros((0, n_dof), dtype=np.float64) 

98 ff_prime_source_function = gsfit_rs.EfitPolynomial(n_dof, regularisations) 

99 elif settings["source_function_ff_prime.json"]["method"] == "tensioned_cubic_b_spline": 

100 ff_prime_tensioned_cubic_b_spline_settings = settings["source_function_ff_prime.json"]["tensioned_cubic_b_spline"] 

101 interior_knots = np.array(ff_prime_tensioned_cubic_b_spline_settings["interior_knots"]) 

102 interval_tensions = np.array(ff_prime_tensioned_cubic_b_spline_settings["interval_tensions"]) 

103 n_dof = len(interior_knots) + 4 

104 if ff_prime_tensioned_cubic_b_spline_settings.get("make_cubic_spline_regularisations", False): 

105 # Build the regularisations from settings (see `tensioned_cubic_splines_regularisations.py`) 

106 regularisations = make_tensioned_cubic_b_spline_regularisations( 

107 interior_knots, 

108 interval_tensions, 

109 ff_prime_tensioned_cubic_b_spline_settings["regularisation_builder"], 

110 ) 

111 else: 

112 regularisations = np.array(ff_prime_tensioned_cubic_b_spline_settings["regularizations"]) 

113 # If `regularisations` is [[]] in the json file, will be interpreted by numpy as having size (1, 0). 

114 # Which would be interpreted as (n_regularisations, n_dof). So it would cause an error 

115 if regularisations.shape == (1, 0): 

116 regularisations = np.zeros((0, n_dof), dtype=np.float64) 

117 ff_prime_source_function = gsfit_rs.TensionedCubicBSpline(regularisations, interior_knots, interval_tensions) 

118 else: 

119 raise ValueError(f"Unknown method for ff_prime source function: {settings['source_function_ff_prime.json']['method']}") 

120 

121 # Grid size and shape 

122 n_r = settings["GSFIT_code_settings.json"]["grid"]["n_r"] 

123 n_z = settings["GSFIT_code_settings.json"]["grid"]["n_z"] 

124 r_min = settings["GSFIT_code_settings.json"]["grid"]["r_min"] 

125 r_max = settings["GSFIT_code_settings.json"]["grid"]["r_max"] 

126 z_min = settings["GSFIT_code_settings.json"]["grid"]["z_min"] 

127 z_max = settings["GSFIT_code_settings.json"]["grid"]["z_max"] 

128 

129 # Normalised poloidal flux grid 

130 n_psi_n = settings["GSFIT_code_settings.json"]["n_psi_n"] 

131 psi_n = np.linspace(0.0, 1.0, n_psi_n).astype(np.float64) 

132 

133 # Initialise the Plasma Rust class 

134 plasma = Plasma( 

135 n_r, 

136 n_z, 

137 r_min, 

138 r_max, 

139 z_min, 

140 z_max, 

141 psi_n, # BUXTON: perhaps better to send in `n_psi_n` 

142 p_prime_source_function, 

143 ff_prime_source_function, 

144 initial_guess_ip, 

145 initial_guess_cur_r, 

146 initial_guess_cur_z, 

147 initial_guess_minor_radius, 

148 initial_guess_elongation, 

149 n_iter_max, 

150 n_iter_min, 

151 n_iter_no_vertical_feedback, 

152 grad_shafranov_deviation_tolerance, 

153 use_anderson_mixing, 

154 anderson_mixing_from_previous_iter, 

155 times_to_reconstruct, 

156 pressure_edge_free, 

157 pressure_edge_regularisation_weight, 

158 ) 

159 

160 return plasma