Coverage for src/swiift/lib/potential.py: 0%

1import itertools 

2 

3import attrs 

4import numpy as np 

5import scipy.sparse as sparse 

6 

7from ..model.model import DiscreteSpectrum, FreeSurfaceWaves, WavesUnderFloe 

8 

9 

10@attrs.define 

11class FilterSpectrum: 

12 gravity: float 

13 spectrum: DiscreteSpectrum 

14 fsw: FreeSurfaceWaves 

15 wuf: WavesUnderFloe 

16 

17 @classmethod 

18 def from_ice_ocean(cls): 

19 pass 

20 

21 

22def scattering(self, floes: Sequence[Floe]) -> list[FloeCoupled]: 

23 """Compute the scattering due to the array of floe. 

24 

25 Each floe in the domain impacts the propagation of the far-field 

26 wave. Every edge causes some reflection and transmission. 

27 This method compute the back-propagating and forward-propagating 

28 waves generated by these edges for a given sequence of `Floe` 

29 and returns the corresponding sequence of `FloeCoupled`. 

30 

31 Paramerers 

32 ---------- 

33 floes : Sequence of Floe 

34 Array of floes causing the scattering. 

35 

36 Returns 

37 ------- 

38 list of FloeCoupled 

39 The floes in the domain with fully-determined radiated fields. 

40 

41 """ 

42 if len(floes) >= 1: 

43 

44 def water_phase(j: int) -> np.ndarray: 

45 assert 1 <= j <= nf - 1 

46 return np.exp( 

47 1j 

48 * self.ocean.wave_numbers[:nvp1] 

49 * (floes[j].left_edge - (floes[j - 1].left_edge + floes[j - 1].length)) 

50 ) 

51 

52 def ice_phase(j: int) -> np.ndarray: 

53 assert 0 <= j <= nf - 1 

54 floe = floes[j] 

55 if not np.isinf(floe.length): 

56 return np.exp(1j * floe.length * floe.ice.wave_numbers[:nvp3]) 

57 else: 

58 # si la longueur du floe est infinie, toutes les 

59 # atténuations associées à des nombres d'ondes avec une 

60 # partie imaginaire non-nulle tendent vers 0 

61 msk = np.imag(floe.ice.wave_numbers[:nvp3]) == 0 

62 phase = np.real(floe.ice.wave_numbers[:nvp3]) % (2 * PI) 

63 return np.exp(1j * phase) * msk 

64 

65 nf = len(floes) 

66 bs = 4 * self.nv + 8 # size of a square block in lhs matrix 

67 nvp3 = self.nv + 3 

68 nvp1 = self.nv + 1 

69 hbs = bs // 2 

70 nfm1 = nf - 1 

71 nfp1 = nf + 1 

72 

73 # Cinq étapes: 

74 # départ avec une matrice identité; 

75 # remplissage de la diagonale-bloc principale, en haut puis 

76 # en bas, le padding à droite (res à gauche) étant ajouté 

77 # avec un produit matriciel; 

78 # construction des 2ndes diagonales-bloc supérieure puis 

79 # inférieure, le padding étant également ajouté via un 

80 # produit matriciel. 

81 

82 # True diagonal 

83 lhs = sparse.eye(nf * bs, format="csr") 

84 # principal block-diagonal, above the diagonal 

85 arr = np.empty((nf, hbs, nvp3), complex) 

86 for i in range(nf): 

87 ice = floes[i].ice 

88 arr[i, :nvp1, :] = -ice.transmission_itow * ice_phase(i) 

89 arr[i, nvp1:, :] = -ice.reflection_i * ice_phase(i) 

90 indices = np.arange(0, nf) 

91 indptr = np.zeros(2 * nf + 1) 

92 indptr[:-1:2] = np.arange(0, nf) 

93 indptr[1::2] = np.arange(1, nfp1) 

94 indptr[-1] = nf 

95 indices1 = np.arange(1, 2 * nf, 2) 

96 indptr1 = np.arange(0, nfp1) 

97 diag = sparse.bsr_matrix( 

98 ( 

99 np.broadcast_to(np.eye(nvp3, hbs, dtype="int8"), (nf, nvp3, hbs)), 

100 indices1, 

101 indptr1, 

102 ), 

103 shape=(nvp3 * nf, bs * nf), 

104 ).tocsr() 

105 diag.eliminate_zeros() 

106 lhs += ( 

107 sparse.bsr_matrix( 

108 (arr, indices, indptr), shape=(nf * bs, nf * nvp3) 

109 ).tocsr() 

110 @ diag 

111 ) 

112 # principal block-diagonal, below the diagonal 

113 for i in range(nf): 

114 ice = floes[i].ice 

115 arr[i, :nvp3, :] = -ice.reflection_i * ice_phase(i) 

116 arr[i, nvp3:, :] = -ice.transmission_itow * ice_phase(i) 

117 indptr[1:] = indptr[0:-1] 

118 indptr[0] = 0 

119 diag = sparse.bsr_matrix( 

120 ( 

121 np.broadcast_to(np.eye(nvp3, hbs, nvp1, dtype="int8"), (nf, nvp3, hbs)), 

122 indices1 - 1, 

123 indptr1, 

124 ), 

125 shape=(nvp3 * nf, bs * nf), 

126 ).tocsr() 

127 diag.eliminate_zeros() 

128 lhs += ( 

129 sparse.bsr_matrix( 

130 (arr, indices, indptr), shape=(nf * bs, nf * nvp3) 

131 ).tocsr() 

132 @ diag 

133 ) 

134 # 2nd block-diagonal, upper 

135 arr = np.empty((nfm1, hbs, nvp1), complex) 

136 for i in range(nfm1): 

137 ice = floes[i].ice 

138 arr[i, :nvp3, :] = -ice.transmission_wtoi * water_phase(i + 1) 

139 arr[i, nvp3:, :] = -ice.reflection_w * water_phase(i + 1) 

140 indices = np.arange(1, nf) 

141 indptr = np.zeros(2 * nf + 1, dtype=int) 

142 indptr[0:-1:2] = np.arange(0, nf) 

143 indptr[1::2] = np.arange(0, nf) 

144 indptr[-1] = nfm1 

145 indices1 = np.arange(2, 2 * nf, 2) 

146 indptr1 = np.zeros(nfp1, dtype="int") 

147 indptr1[1:] = np.arange(0, nf) 

148 diag = sparse.bsr_matrix( 

149 ( 

150 np.broadcast_to(np.eye(nvp1, hbs, dtype="int8"), (nfm1, nvp1, hbs)), 

151 indices1, 

152 indptr1, 

153 ), 

154 shape=(nvp1 * nf, bs * nf), 

155 ).tocsr() 

156 diag.eliminate_zeros() 

157 lhs += ( 

158 sparse.bsr_matrix( 

159 (arr, indices, indptr), shape=(nf * bs, nf * nvp1) 

160 ).tocsr() 

161 @ diag 

162 ) 

163 # 2nd block-diagonal, lower 

164 for i in range(nfm1): 

165 ice = floes[i + 1].ice 

166 arr[i, :nvp1, :] = -ice.reflection_w * water_phase(i + 1) 

167 arr[i, nvp1:, :] = -ice.transmission_wtoi * water_phase(i + 1) 

168 indptr[1:] = indptr[0:-1] 

169 indptr[0] = 0 

170 indptr1[0:-1] = indptr1[1:] 

171 indptr1[-1] = nfm1 

172 diag = sparse.bsr_matrix( 

173 ( 

174 np.broadcast_to( 

175 np.eye(nvp1, hbs, nvp3, dtype="int8"), (nfm1, nvp1, hbs) 

176 ), 

177 indices1 - 1, 

178 indptr1, 

179 ), 

180 shape=(nvp1 * nf, bs * nf), 

181 ).tocsr() 

182 diag.eliminate_zeros() 

183 lhs += ( 

184 sparse.bsr_matrix( 

185 (arr, indices - 1, indptr), shape=(nf * bs, nf * nvp1) 

186 ).tocsr() 

187 @ diag 

188 ) 

189 

190 # forcing; it is assumed that nothing comes from the right 

191 arr = np.empty((hbs, 1), dtype=complex) 

192 arr[:nvp1, 0] = floes[0].ice.reflection_w[:, 0] 

193 arr[nvp1:hbs, 0] = floes[0].ice.transmission_wtoi[:, 0] 

194 arr *= self.wave.potential_amplitude 

195 rhs = sparse.csc_matrix(arr) 

196 

197 # si le dernier floe est semi-infini, cas spécial : 

198 # le dernier edge est une transition eau vers glace, 

199 # un demi-bloc de moins est nécessaire. 

200 if floes[-1].length == np.inf: 

201 lhs = lhs[: nfm1 * bs + hbs, : nfm1 * bs + hbs] 

202 rhs.resize(((nfm1 * bs + hbs, 1))) 

203 else: 

204 rhs.resize(((nf * bs, 1))) 

205 

206 # résolution du système 

207 # spsolve options from benchmark on laptop 

208 sol = sparse.linalg.spsolve(lhs, rhs, "NATURAL", False) 

209 

210 ww_neg = (sol[i * bs : i * bs + nvp1] for i in range(nf)) 

211 iw_pos = (sol[i * bs + nvp1 : i * bs + hbs] for i in range(nf)) 

212 if floes[-1].length == np.inf: 

213 # No right_edge, so it doesn't radiate anything 

214 l_iw_neg = [None] 

215 l_ww_pos = [None] 

216 else: 

217 # Contribution of the right edge of the rightmost floe 

218 l_iw_neg = [sol[nfm1 * bs + hbs : nfm1 * bs + hbs + nvp3]] 

219 l_ww_pos = [sol[nfm1 * bs + hbs + nvp3 : nf * bs]] 

220 iw_neg = (sol[i * bs + hbs : i * bs + hbs + nvp3] for i in range(nfm1)) 

221 ww_pos = (sol[i * bs + hbs + nvp3 : (i + 1) * bs] for i in range(nfm1)) 

222 

223 return [ 

224 FloeCoupled(floe, (iw_n, iw_p, ww_n, ww_p)) 

225 for (floe, iw_n, iw_p, ww_n, ww_p) in zip( 

226 floes, 

227 itertools.chain(iw_neg, l_iw_neg), 

228 iw_pos, 

229 ww_neg, 

230 itertools.chain(ww_pos, l_ww_pos), 

231 ) 

232 ]