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 

4 

5 

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

7 """ 

8 Calculate the Greens table with the boundary points. 

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

10 

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

12 

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

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

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

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

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

18 

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

20 """ 

21 

22 r = plasma.get_array1(["grid", "r"]) 

23 z = plasma.get_array1(["grid", "z"]) 

24 r_min = np.min(r) 

25 r_max = np.max(r) 

26 z_min = np.min(z) 

27 z_max = np.max(z) 

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

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

30 ones_len_r_inner = np.ones_like(r_inner) 

31 ones_len_z_inner = np.ones_like(z_inner) 

32 

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

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

35 

36 r_ltrb = np.concatenate( 

37 ( 

38 [r_min], # (bottom, left) 

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

40 [r_min], # (top, left) 

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

42 [r_max], # (top, right) 

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

44 [r_max], # (bottom, right) 

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

46 ) 

47 ) 

48 z_ltrb = np.concatenate( 

49 ( 

50 [z_min], # (bottom, left) 

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

52 [z_max], # (top, left) 

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

54 [z_max], # (top, right) 

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

56 [z_min], # (bottom, right) 

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

58 ) 

59 ) 

60 

61 n_ltrb = len(r_ltrb) 

62 

63 d_r_vec = np.ones(n_ltrb) * d_r 

64 d_z_vec = np.ones(n_ltrb) * d_z 

65 

66 # Calculate the inductance matrix between boundary points 

67 g_ltrb = gsfit_rs.greens_py( 

68 r_ltrb, 

69 z_ltrb, 

70 r_ltrb, 

71 z_ltrb, 

72 d_r_vec, # Needed for self-inductance 

73 d_z_vec, # Needed for self-inductance 

74 ) 

75 

76 # Check diagonal values for self-inductance are g_ltrb[i, i] = self_inductance_rectangle_cross_section(r_ltrb[i], d_r, d_z) 

77 for i in range(n_ltrb): 

78 expected_self_inductance = self_inductance_rectangle_cross_section(r_ltrb[i], d_r_vec[i], d_z_vec[i]) 

79 if not np.isclose(g_ltrb[i, i], expected_self_inductance): 

80 raise ValueError(f"Diagonal value at index {i} does not match expected self-inductance: {g_ltrb[i, i]} != {expected_self_inductance}") 

81 

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

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

84 # Note the corners are at the following indicdes: 

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

86 tl_corner_idx = len(z) - 1 

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

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

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

90 # Left side 

91 for i in range(tl_corner_idx): 

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

93 # Top side 

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

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

96 # Right side 

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

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

99 # Bottom side 

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

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

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

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

104 

105 # Divide through by 2pi 

106 g_ltrb /= 2 * np.pi 

107 

108 return g_ltrb.flatten() 

109 

110 

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

112 """ 

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

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

115 

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

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

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

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

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

121 """ 

122 

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

124 return ( 

125 mu_0 

126 * r 

127 * ( 

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

129 - 0.5 

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

131 ) 

132 )