Coverage for python/gsfit/database_writers/rtgsfit_mdsplus/greens_with_boundary_points.py: 0%

1import gsfit_rs 

2import numpy as np 

3import numpy.typing as npt 

4from gsfit_rs.imas import equilibrium_paths as ep 

5 

6 

7def greens_with_boundary_points(plasma: gsfit_rs.Plasma) -> npt.NDArray[np.float64]: 

8 """ 

9 Calculate the Greens table with the boundary points. 

10 In this context "boundary" means the 2D (R, Z) grid, not the plasma boundary. 

11 

12 :param plasma: Plasma object containing the grid information. 

13 

14 Boundary points are ordered clockwise from the bottom left, i.e. 

15 * (bottom, left) = (r_min, z_min) 

16 * (top, left) = (r_min, z_max) 

17 * (top, right) = (r_max, z_max) 

18 * (bottom, right) = (r_max, z_min) 

19 

20 Note: need to be careful not to double count the corner points 

21 """ 

22 

23 # `profiles_2d(0)` because GSFit solves on a single rectangular (R, Z) grid 

24 r = plasma.get(ep.time_slice[0].profiles_2d[0].grid.dim1) 

25 z = plasma.get(ep.time_slice[0].profiles_2d[0].grid.dim2) 

26 r_min = np.min(r) 

27 r_max = np.max(r) 

28 z_min = np.min(z) 

29 z_max = np.max(z) 

30 r_inner = r[1:-1] # Exclude the first and last points to avoid double counting corners 

31 z_inner = z[1:-1] # Exclude the first and last points to avoid double counting corners 

32 ones_len_r_inner = np.ones_like(r_inner) 

33 ones_len_z_inner = np.ones_like(z_inner) 

34 

35 d_r = np.mean(np.diff(r)).astype(np.float64) 

36 d_z = np.mean(np.diff(z)).astype(np.float64) 

37 

38 r_ltrb = np.concatenate( 

39 ( 

40 [r_min], # (bottom, left) 

41 r_min * ones_len_z_inner, # traverse (bottom, left)-delta to (top, left)-delta 

42 [r_min], # (top, left) 

43 r_inner, # traverse (top, left)-delta to (top, right)-delta 

44 [r_max], # (top, right) 

45 r_max * ones_len_z_inner, # traverse (top, right)-delta to (bottom, right)-delta 

46 [r_max], # (bottom, right) 

47 np.flip(r_inner), # traverse (bottom, right)-delta to (bottom, left)-delta 

48 ) 

49 ) 

50 z_ltrb = np.concatenate( 

51 ( 

52 [z_min], # (bottom, left) 

53 z_inner, # traverse (bottom, left)-delta to (top, left)-delta 

54 [z_max], # (top, left) 

55 z_max * ones_len_r_inner, # traverse (top, left)-delta to (top, right)-delta 

56 [z_max], # (top, right) 

57 np.flip(z_inner), # traverse (top, right)-delta to (bottom, right)-delta 

58 [z_min], # (bottom, right) 

59 z_min * ones_len_r_inner, # traverse (bottom, right)-delta to (bottom, left)-delta 

60 ) 

61 ) 

62 

63 n_ltrb = len(r_ltrb) 

64 

65 d_r_vec = np.ones(n_ltrb) * d_r 

66 d_z_vec = np.ones(n_ltrb) * d_z 

67 

68 # Calculate the inductance matrix between boundary points 

69 g_ltrb = gsfit_rs.greens_py( 

70 r_ltrb, 

71 z_ltrb, 

72 r_ltrb, 

73 z_ltrb, 

74 d_r_vec, # Needed for self-inductance 

75 d_z_vec, # Needed for self-inductance 

76 ) 

77 

78 # The diagonal self-inductance values returned by greens_py are evaluated at the centre of the 

79 # rectangular cross-section. Here we instead need the average over the cross-section, since g_ltrb 

80 # is used to compute a line integral for which the cross-section average is more accurate. We 

81 # therefore overwrite the diagonals below rather than validating them against greens_py. 

82 

83 # Need to overwrite the diagonals of the matrix with the self-inductance for a surface current on 

84 # the computational boundary. Excluding the corners, which are not used by RTGSFIT 

85 # Note the corners are at the following indicdes: 

86 # 0, len(z) - 1, len(z) + len(r) - 2, 2 * len(z) + len(r) - 3, 2 * len(z) + 2 * len(r) - 4 

87 tl_corner_idx = len(z) - 1 

88 tr_corner_idx = len(z) + len(r) - 2 

89 br_corner_idx = 2 * len(z) + len(r) - 3 

90 bl_corner_idx = 2 * len(z) + 2 * len(r) - 4 

91 # Left side 

92 for i in range(tl_corner_idx): 

93 g_ltrb[i, i] = self_inductance_rectangle_cross_section(r_ltrb[i], 0, d_z) 

94 # Top side 

95 for i in range(tl_corner_idx + 1, tr_corner_idx): 

96 g_ltrb[i, i] = self_inductance_rectangle_cross_section(r_ltrb[i], d_r, 0) 

97 # Right side 

98 for i in range(tr_corner_idx + 1, br_corner_idx): 

99 g_ltrb[i, i] = self_inductance_rectangle_cross_section(r_ltrb[i], 0, d_z) 

100 # Bottom side 

101 for i in range(br_corner_idx + 1, bl_corner_idx): 

102 g_ltrb[i, i] = self_inductance_rectangle_cross_section(r_ltrb[i], d_r, 0) 

103 # Note: g_ltrb.shape = (n_ltrb, n_ltrb) = (2 * n_r + 2 * n_z - 4, 2 * n_r + 2 * n_z - 4) 

104 # The -4 is so that we don't double count corners 

105 

106 # Divide through by 2pi 

107 g_ltrb /= 2 * np.pi 

108 

109 return g_ltrb.flatten() 

110 

111 

112def self_inductance_rectangle_cross_section(r: float, delta_r: float, delta_z: float) -> float: 

113 """ 

114 Calculate the self-inductance of an axisymmetric wire with a rectangular cross-section 

115 of height delta_z and width delta_r, centered at (r, z). 

116 

117 :param r: The R coordinate of the center of the rectangle. 

118 :param z: The Z coordinate of the center of the rectangle. 

119 :param delta_r: The width of the rectangle in the R direction. 

120 :param delta_z: The height of the rectangle in the Z direction. 

121 :return: The self-inductance of the rectangle in H (Henries). 

122 """ 

123 

124 mu_0 = 4e-7 * np.pi # Vacuum permeability in H/m 

125 return ( 

126 mu_0 

127 * r 

128 * ( 

129 (1 + 2 * (delta_z / (8 * r)) ** 2 + 2 / 3 * (delta_r / (8 * r)) ** 2) * float(np.log(8 * r / (delta_r + delta_z))) 

130 - 0.5 

131 + 0.5 * (delta_z / (8 * r)) ** 2 

132 ) 

133 )