Coverage for src/swiift/lib/physics.py: 35%
233 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-09-11 16:23 +0200
1from __future__ import annotations
3import functools
4import typing
6import attrs
7import numpy as np
9from . import numerical
10from ._ph_utils import _wavefield
11from .constants import PI_D4, SQR2
13if typing.TYPE_CHECKING:
14 from ..model import model
17def _package_wuf(
18 wuf: model.WavesUnderFloe, growth_params: tuple | None
19) -> tuple[tuple, tuple, tuple | None]:
20 floe_params = wuf.wui.ice._red_elastic_number, wuf.length
21 wave_params = wuf.edge_amplitudes, wuf.wui._c_wavenumbers
22 if growth_params is not None:
23 growth_params = growth_params[0] - wuf.left_edge, growth_params[1]
24 return floe_params, wave_params, growth_params
27def _dis_par_amps(red_num: float, wave_params: tuple[np.ndarray, np.ndarray]):
28 """Complex amplitudes of individual particular solutions"""
29 c_amplitudes, c_wavenumbers = wave_params
30 return c_amplitudes / (1 + 0.25 * (c_wavenumbers / red_num) ** 4)
33def _dis_hom_coefs(
34 floe_params: tuple[float, float], wave_params: tuple[np.ndarray, np.ndarray]
35) -> np.ndarray:
36 """Coefficients of the four orthogonal homogeneous solutions"""
37 return _dis_hom_mat(*floe_params) @ _dis_hom_rhs(floe_params, wave_params)
40def _dis_hom_mat(red_num: float, length: float):
41 """Linear application to determine, from the BCs, the coefficients of
42 the four independent solutions to the homo ODE"""
43 adim = red_num * length
44 adim2 = 2 * adim
45 denom = (
46 -2
47 * red_num**2
48 * (np.expm1(-adim2) ** 2 + 2 * np.exp(-adim2) * (np.cos(adim2) - 1))
49 )
50 mat = np.array(
51 [
52 [
53 SQR2 * np.exp(-adim) * np.sin(adim2 + PI_D4) - np.exp(-3 * adim),
54 -SQR2 * np.cos(adim + PI_D4)
55 - 3 * np.exp(-adim2) * np.sin(adim)
56 + np.exp(-adim2) * np.cos(adim),
57 np.exp(-adim) * np.sin(adim2) - np.exp(-adim) + np.exp(-3 * adim),
58 np.cos(adim)
59 - 2 * np.exp(-adim2) * np.sin(adim)
60 - np.exp(-adim2) * np.cos(adim),
61 ],
62 [
63 -SQR2 * np.exp(-adim) * np.cos(adim2 + PI_D4)
64 + 2 * np.exp(-adim)
65 - np.exp(-3 * adim),
66 -SQR2 * np.sin(adim + PI_D4)
67 + SQR2 * np.exp(-adim2) * np.cos(adim + PI_D4),
68 -np.exp(-adim) * np.cos(adim2) + np.exp(-adim),
69 np.sin(adim) - np.exp(-adim2) * np.sin(adim),
70 ],
71 [
72 -1 + SQR2 * np.exp(-adim2) * np.cos(adim2 + PI_D4),
73 (3 * np.sin(adim) + np.cos(adim)) * np.exp(-adim)
74 - SQR2 * np.exp(-3 * adim) * np.sin(adim + PI_D4),
75 (np.sin(adim2) + 1) * np.exp(-adim2) - 1,
76 (-2 * np.sin(adim) + np.cos(adim)) * np.exp(-adim)
77 - np.exp(-3 * adim) * np.cos(adim),
78 ],
79 [
80 (SQR2 * np.sin(adim2 + PI_D4) - 2) * np.exp(-adim2) + 1,
81 -SQR2 * np.exp(-adim) * np.sin(adim + PI_D4)
82 + SQR2 * np.exp(-3 * adim) * np.cos(adim + PI_D4),
83 (-np.cos(adim2) + 1) * np.exp(-adim2),
84 np.exp(-adim) * np.sin(adim) - np.exp(-3 * adim) * np.sin(adim),
85 ],
86 ]
87 )
88 mat[:, 2:] /= red_num
90 return mat / denom
93def _dis_hom_rhs(
94 floe_params: tuple[float, float], wave_params: tuple[np.ndarray, np.ndarray]
95):
96 """Vector onto which apply the matrix, to extract the coefficients"""
97 red_num, length = floe_params
98 _, c_wavenumbers = wave_params
99 exp_arg = 1j * c_wavenumbers * length
101 r1 = c_wavenumbers**2 * _dis_par_amps(red_num, wave_params)
102 r2 = np.imag(r1 @ np.exp(exp_arg))
103 r3 = 1j * c_wavenumbers * r1
104 r4 = np.imag(r3 @ np.exp(exp_arg))
105 r1 = np.imag(r1).sum()
106 r3 = np.imag(r3).sum()
108 return np.array((r1, r2, r3, r4))
111# Can be used to decorate methods to make them return a scalar instead of
112# a 1d-array of length 1, if the argument is a scalar OR a 0d array (that is, a
113# scalar). If the argument is a tuple or a list or a 1d array of length 1,
114# a 1d-array of length 1 is returned.
115# This emulate the behaviour of standard numpy ufuncs.
116def _demote_to_scalar(f):
117 @functools.wraps(f)
118 def wrapper(*args):
119 instance, x = args 1bc
120 x = np.asarray(x) 1bc
121 res = f(instance, x) 1bc
122 if len(res) == 1 and x.ndim == 0: 1bc
123 return res.item() 1bc
124 return res 1b
126 return wrapper
129@attrs.define
130class FluidSurfaceHandler:
131 wave_params: tuple[np.ndarray, np.ndarray]
132 growth_params: tuple[np.ndarray, float] | None = None
134 @classmethod
135 def from_wuf(cls, wuf: model.WavesUnderFloe, growth_params=None):
136 return cls(*(_package_wuf(wuf, growth_params)[1:]))
138 @classmethod
139 def from_domain(cls, domain: model.Domain, growth_params: tuple | None = None):
140 return cls(
141 (
142 domain.spectrum.amplitudes * np.exp(1j * domain.spectrum.phases),
143 domain.fsw.wavenumbers,
144 ),
145 growth_params,
146 )
148 @_demote_to_scalar
149 def compute(self, x, /):
150 return numerical.free_surface(x, self.wave_params, self.growth_params)
153@attrs.define
154class DisplacementHandler:
155 floe_params: tuple[float, float]
156 wave_params: tuple[np.ndarray, np.ndarray]
157 growth_params: tuple[np.ndarray, float] | None = None
159 @classmethod
160 def from_wuf(cls, wuf: model.WavesUnderFloe, growth_params=None):
161 return cls(*_package_wuf(wuf, growth_params))
163 def _dis_hom(self, x: np.ndarray):
164 """Homogeneous solution to the displacement ODE"""
165 red_num, length = self.floe_params
166 arr = red_num * x
167 cosx, sinx = np.cos(arr), np.sin(arr)
168 expx = np.exp(-red_num * (length - x))
169 exmx = np.exp(-arr)
170 return np.vstack(
171 ([expx * cosx, expx * sinx, exmx * cosx, exmx * sinx])
172 ).T @ _dis_hom_coefs(self.floe_params, self.wave_params)
174 def _dis_par(self, x: np.ndarray):
175 """Sum of the particular solutions to the displacement ODE"""
176 return _wavefield(
177 x, _dis_par_amps(self.floe_params[0], self.wave_params), self.wave_params[1]
178 )
180 @_demote_to_scalar
181 def _dis(self, x, /):
182 return self._dis_hom(x) + self._dis_par(x)
184 def compute(
185 self,
186 x,
187 an_sol: bool | None = None,
188 num_params: dict | None = None,
189 ):
190 if numerical._use_an_sol(
191 an_sol, self.floe_params[1], self.growth_params, linear_curvature=None
192 ):
193 return self._dis(np.asarray(x))
194 return numerical.displacement(
195 x, self.floe_params, self.wave_params, self.growth_params, num_params
196 )
199@attrs.define
200class CurvatureHandler:
201 floe_params: tuple[float, float]
202 wave_params: tuple[np.ndarray, np.ndarray]
203 growth_params: tuple[np.ndarray, float] | None = None
205 @classmethod
206 def from_wuf(cls, wuf: model.WavesUnderFloe, growth_params=None):
207 return cls(*_package_wuf(wuf, growth_params))
209 def _cur_wavefield(self, x):
210 """Second derivative of the interface"""
211 red_num = self.floe_params[0]
212 _, c_wavenumbers = self.wave_params
214 return -np.imag(
215 (_dis_par_amps(red_num, self.wave_params) * c_wavenumbers**2)
216 @ np.exp(1j * c_wavenumbers[:, None] * x)
217 )
219 def _cur_hom(self, x):
220 """Second derivative of the homogeneous part of the displacement"""
221 red_num, length = self.floe_params
222 arr = red_num * x
223 cosx, sinx = np.cos(arr), np.sin(arr)
224 expx = np.exp(-red_num * (length - x))
225 exmx = np.exp(-arr)
226 return (
227 2
228 * red_num**2
229 * np.vstack(([-expx * sinx, expx * cosx, exmx * sinx, -exmx * cosx])).T
230 @ _dis_hom_coefs(self.floe_params, self.wave_params)
231 )
233 def _cur_par(self, x):
234 """Second derivative of the particular part of the displacement"""
235 return self._cur_wavefield(x)
237 @_demote_to_scalar
238 def _cur(self, x, /):
239 return self._cur_hom(x) + self._cur_par(x)
241 def compute(
242 self,
243 x,
244 an_sol: bool | None = None,
245 num_params: dict | None = None,
246 is_linear: bool = True,
247 ):
248 """Curvature of the floe.
250 A linear approximation is given by the second derivative of the
251 displacement. The exact definition involves the first derivative of the
252 displacement.
254 Parameters
255 ----------
256 x : float or 1d array_like of float
257 Along-floe coordinates at which to compute the curvature, in m
258 an_sol : bool
259 Whether to use an analytical formulation of the curvature
260 num_params : dict
261 Optional arguments passed to the numerical solver if `not an_sol`
262 is_linear : bool
263 Whether to use the linear approximation; ignored if `an_sol`
265 Returns
266 -------
267 float or 1d array_like of float
269 """
270 if numerical._use_an_sol(
271 an_sol, self.floe_params[1], self.growth_params, is_linear
272 ):
273 return self._cur(x)
274 return numerical.curvature(
275 x,
276 self.floe_params,
277 self.wave_params,
278 self.growth_params,
279 num_params,
280 is_linear,
281 )
284@attrs.define
285class StrainHandler:
286 curv_handler: CurvatureHandler
287 thickness: float
289 @classmethod
290 def from_wuf(cls, wuf: model.WavesUnderFloe, growth_params=None):
291 return cls(CurvatureHandler.from_wuf(wuf, growth_params), wuf.wui.ice.thickness)
293 def compute(self, x, an_sol: bool | None = None, num_params: dict | None = None):
294 return -self.thickness / 2 * self.curv_handler.compute(x, an_sol, num_params)
297@attrs.define
298class EnergyHandler:
299 floe_params: tuple[float, float]
300 wave_params: tuple[np.ndarray, np.ndarray]
301 growth_params: tuple[np.ndarray, float] | None = None
302 factor: float = 1
304 @classmethod
305 def from_wuf(cls, wuf: model.WavesUnderFloe, growth_params=None):
306 factor = cls._compute_factor(wuf)
307 return cls(*_package_wuf(wuf, growth_params), factor)
309 @staticmethod
310 def _compute_factor(wuf: model.WavesUnderFloe):
311 # NOTE: this could also be a WUF property?
312 return wuf.wui.ice.flex_rigidity / (2 * wuf.wui.ice.thickness)
314 def _egy_hom(self):
315 """Energy from the homogen term of the displacement ODE"""
316 red_num, length = self.floe_params
317 adim = red_num * length
318 adim2 = 2 * adim
319 c_1, c_2, c_3, c_4 = _dis_hom_coefs(self.floe_params, self.wave_params)
321 return red_num**3 * (
322 +(c_1**2) * (-SQR2 * np.sin(adim2 + PI_D4) + 2 - np.exp(-adim2)) / 2
323 + c_2**2 * (SQR2 * np.sin(adim2 + PI_D4) + 2 - 3 * np.exp(-adim2)) / 2
324 + c_3**2 * ((SQR2 * np.cos(adim2 + PI_D4) - 2) * np.exp(-adim2) + 1) / 2
325 + c_4**2 * (3 - (SQR2 * np.cos(adim2 + PI_D4) + 2) * np.exp(-adim2)) / 2
326 + c_1 * c_2 * (SQR2 * np.cos(adim2 + PI_D4) - np.exp(-adim2))
327 - 2
328 * red_num
329 * c_1
330 * c_3
331 * (2 * length - np.sin(adim2) / red_num)
332 * np.exp(-adim)
333 + 4 * c_1 * c_4 * np.exp(-adim) * np.sin(adim) ** 2
334 + 4 * c_2 * c_3 * np.exp(-adim) * np.sin(adim) ** 2
335 - 2
336 * red_num
337 * c_2
338 * c_4
339 * (2 * length + np.sin(adim2) / red_num)
340 * np.exp(-adim)
341 + c_3 * c_4 * (-1 + SQR2 * np.exp(-adim2) * np.sin(adim2 + PI_D4))
342 )
344 def _egy_par_vals(self):
345 red_num = self.floe_params[0]
346 comp_amps = _dis_par_amps(red_num, self.wave_params)
347 _, c_wavenumbers = self.wave_params
349 return comp_amps * (1j * c_wavenumbers) ** 2
351 def _egy_par_pow2(self):
352 """Energy contribution from individual forcings"""
353 red_num, length = self.floe_params
354 _, c_wavenumbers = self.wave_params
355 wavenumbers = np.real(c_wavenumbers)
356 attenuations = np.imag(c_wavenumbers)
358 comp_curvs = self._egy_par_vals()
359 wn_moduli, curv_moduli = map(np.abs, (c_wavenumbers, comp_curvs))
360 wn_phases, curv_phases = map(np.angle, (c_wavenumbers, comp_curvs))
362 # Because of a division by `attenuations`, this term must me handled
363 # with some extra care
364 attenuation_mask = attenuations != 0
365 non_oscillatory_term = np.full(attenuations.shape, np.nan)
366 non_oscillatory_term[attenuation_mask] = (
367 -np.expm1(-2 * attenuations[attenuation_mask] * length)
368 / attenuations[attenuation_mask]
369 )
370 # NOTE: corresponds to the limit of the previous term, when attenuation -> 0
371 non_oscillatory_term[~attenuation_mask] = 2 * length
373 return (
374 curv_moduli**2
375 @ (
376 non_oscillatory_term
377 + (
378 np.sin(2 * curv_phases - wn_phases)
379 - np.sin(2 * (wavenumbers * length + curv_phases) - wn_phases)
380 * np.exp(-2 * attenuations * length)
381 )
382 / wn_moduli
383 )
384 ) / 4
386 def _egy_par_m(self):
387 """Energy contribution from forcing interactions"""
388 red_num, length = self.floe_params
389 _, c_wavenumbers = self.wave_params
390 wavenumbers = np.real(c_wavenumbers)
391 attenuations = np.imag(c_wavenumbers)
392 comp_curvs = self._egy_par_vals()
394 # Binomial coefficients, much quicker than itertools
395 idx1, idx2 = np.triu_indices(c_wavenumbers.size, 1)
397 mean_attenuations = attenuations[idx1] + attenuations[idx2]
398 comp_wns = (
399 wavenumbers[idx1] - wavenumbers[idx2],
400 wavenumbers[idx1] + wavenumbers[idx2],
401 )
402 comp_wns += 1j * mean_attenuations
404 curv_moduli = np.abs(comp_curvs)
405 _curv_phases = np.angle(comp_curvs)
406 curv_phases = (
407 _curv_phases[idx1] - _curv_phases[idx2],
408 _curv_phases[idx1] + _curv_phases[idx2],
409 )
411 def _f(comp_wns, curv_phases):
412 wn_moduli = np.abs(comp_wns)
413 wn_phases = np.angle(comp_wns)
414 return (
415 np.sin(curv_phases - wn_phases)
416 - (
417 np.exp(-np.imag(comp_wns) * length)
418 * np.sin(np.real(comp_wns) * length + curv_phases - wn_phases)
419 )
420 ) / wn_moduli
422 return (
423 curv_moduli[idx1]
424 * curv_moduli[idx2]
425 @ (_f(comp_wns[1], curv_phases[1]) - _f(comp_wns[0], curv_phases[0]))
426 )
428 def _egy_par(self):
429 """Energy from the forcing term of the displacement ODE"""
430 return self._egy_par_pow2() + self._egy_par_m()
432 def _egy_m(self):
433 red_num, length = self.floe_params
434 _, c_wavenumbers = self.wave_params
435 adim = red_num * length
436 wavenumbers = np.real(c_wavenumbers)
437 attenuations = np.imag(c_wavenumbers)
439 wn_abs, wn_phases = (_f(c_wavenumbers) for _f in (np.abs, np.angle))
440 comp_curvs = self._egy_par_vals()
441 curv_moduli, curv_phases = np.abs(comp_curvs), np.angle(comp_curvs)
442 num_add = red_num + wavenumbers
443 num_sub = red_num - wavenumbers
444 att_add = red_num + attenuations
445 att_sub = red_num - attenuations
446 q_pp = att_add**2 + num_add**2
447 q_mp = att_sub**2 + num_add**2
448 q_pm = att_add**2 + num_sub**2
449 q_mm = att_sub**2 + num_sub**2
451 phase_deltas = wn_phases - curv_phases
452 sin_phase_deltas = np.sin(phase_deltas)
453 cos_phase_deltas = np.cos(phase_deltas)
454 phase_quads = curv_phases + PI_D4
455 sin_phase_quads = np.sin(phase_quads)
456 cos_phase_quads = np.cos(phase_quads)
457 energ_1_K = np.array(
458 (
459 sin_phase_deltas * (1 / q_mp - 1 / q_mm),
460 cos_phase_deltas * (1 / q_mp + 1 / q_mm),
461 -sin_phase_deltas * (1 / q_pp - 1 / q_pm),
462 -cos_phase_deltas * (1 / q_pp + 1 / q_pm),
463 )
464 )
465 energ_1_b = np.array(
466 (
467 -(sin_phase_quads / q_mp - cos_phase_quads / q_mm),
468 (cos_phase_quads / q_mp - sin_phase_quads / q_mm),
469 -(cos_phase_quads / q_pp - sin_phase_quads / q_pm),
470 -(sin_phase_quads / q_pp - cos_phase_quads / q_pm),
471 )
472 )
473 energ_1 = wn_abs * energ_1_K + SQR2 * red_num * energ_1_b
474 energ_1[:2, :] *= np.exp(-adim)
476 arg_add = num_add * length
477 arg_add_phase = arg_add - phase_deltas
478 sin_arg_ap = np.sin(arg_add_phase)
479 cos_arg_ap = np.cos(arg_add_phase)
480 arg_sub = num_sub * length
481 arg_sub_phase = arg_sub + phase_deltas
482 sin_arg_sp = np.sin(arg_sub_phase)
483 cos_arg_sp = np.cos(arg_sub_phase)
484 arg_add_quad = arg_add + phase_quads
485 sin_arg_aq = np.sin(arg_add_quad)
486 cos_arg_aq = np.cos(arg_add_quad)
487 arg_sub_quad = arg_sub - curv_phases + PI_D4
488 sin_arg_sq = np.sin(arg_sub_quad)
489 cos_arg_sq = np.cos(arg_sub_quad)
490 energ_exp_K = np.array(
491 (
492 sin_arg_ap / q_mp + sin_arg_sp / q_mm,
493 -(cos_arg_ap / q_mp + cos_arg_sp / q_mm),
494 -(sin_arg_ap / q_pp + sin_arg_sp / q_pm),
495 cos_arg_ap / q_pp + cos_arg_sp / q_pm,
496 )
497 )
498 energ_exp_b = np.array(
499 (
500 sin_arg_aq / q_mp - sin_arg_sq / q_mm,
501 -cos_arg_aq / q_mp + cos_arg_sq / q_mm,
502 cos_arg_aq / q_pp - cos_arg_sq / q_pm,
503 sin_arg_aq / q_pp - sin_arg_sq / q_pm,
504 )
505 )
506 energ_exp = wn_abs * energ_exp_K + SQR2 * red_num * energ_exp_b
507 energ_exp[2:, :] *= np.exp(-adim)
508 energ_exp *= np.exp(-attenuations * length)
510 energ = (
511 (energ_1 + energ_exp)
512 @ curv_moduli
513 @ _dis_hom_coefs(self.floe_params, self.wave_params)
514 )
515 return red_num**2 * energ
517 def compute(
518 self,
519 an_sol: bool | None = None,
520 num_params: dict | None = None,
521 integration_method: str | None = None,
522 linear_curvature: bool | None = None,
523 **kwargs,
524 ) -> float:
525 do_use_an_sol = numerical._use_an_sol(
526 an_sol, self.floe_params[1], self.growth_params, linear_curvature
527 )
528 if do_use_an_sol:
529 unit_energy = self._egy_hom() + 2 * self._egy_m() + self._egy_par()
530 else:
531 linear_curvature = True if linear_curvature is None else False
532 # TODO: documenter le comportement par défault.
533 # Inclure un comportement différent en mono (systématique) et poly
534 # (stochastic), avec un warning si l'utilisation d'une méthode
535 # systématique est utilisée en poly, car ça pourrait prendre un
536 # temps infini.
537 unit_energy = numerical.unit_energy(
538 self.floe_params,
539 self.wave_params,
540 self.growth_params,
541 num_params,
542 integration_method,
543 linear_curvature,
544 **kwargs,
545 )
547 return self.factor * unit_energy