Coverage for src/flexfrac1d/model/model.py: 92%
321 statements
« prev ^ index » next coverage.py v7.4.1, created at 2024-08-30 14:00 +0200
1from __future__ import annotations
3from collections.abc import Sequence
4import functools
5import itertools
6from numbers import Real
7import operator
8import typing
10import attrs
11import numpy as np
12from sortedcontainers import SortedList
14from ..lib import att, dr, physics as ph
15from ..lib.constants import PI_2, SQR2
16from ..lib.graphics import plot_displacement
18if typing.TYPE_CHECKING:
19 # Guard against circular imports
20 from . import frac_handlers as fh
23@attrs.define(frozen=True)
24class Wave:
25 """A monochromatic wave.
27 Parameters
28 ----------
29 amplitude : float
30 Wave amplitude in m
31 period : float
32 Wave period in s
33 phase :
34 Wave phase in rad
36 Attributes
37 ----------
38 frequency : float
39 Wave frequency in Hz
40 angular_frequency: float
41 Wave angular frequency in rad s^-1
43 """
45 amplitude: float
46 period: float
47 phase: float = attrs.field(default=0, converter=lambda raw: raw % PI_2)
49 @classmethod
50 def from_frequency(cls, amplitude, frequency, phase=0):
51 """Build a wave from frequency instead of period.
53 Parameters
54 ----------
55 amplitude : float
56 The wave amplitude in m
57 frequency : float
58 The wave frequency in Hz
59 phase :
60 The wave phase in rad
62 Returns
63 -------
64 Wave
66 """
67 return cls(amplitude, 1 / frequency, phase)
69 @functools.cached_property
70 def frequency(self) -> float:
71 return 1 / self.period
73 @functools.cached_property
74 def angular_frequency(self) -> float:
75 return PI_2 / self.period
77 @functools.cached_property
78 def _angular_frequency_pow2(self) -> float:
79 """Return the squared angular frequency.
81 This is a convenience method, this quantity appearing frequently in
82 computations.
84 Returns
85 -------
86 float
89 """
90 return self.angular_frequency**2
93@attrs.define(frozen=True)
94class Ocean:
95 """The fluid bearing ice floes.
97 This class encapsulates the properties of an incompressible ocean of
98 constant depth and given density.
100 Parameters
101 ----------
102 depth : float
103 Ocean depth in m
104 density : float
105 Ocean density in kg m^-3
107 """
109 depth: float = np.inf
110 density: float = 1025
113@attrs.define(frozen=True, eq=False)
114@functools.total_ordering
115class _Subdomain:
116 """A segment localised in space.
118 Parameters
119 ----------
120 left_edge : float
121 Coordinate of the left edge of the domain in m
122 length : float
123 Length of the domain in m
125 Attributes
126 ----------
127 right_edge : float
128 Coordinate of the right edge of the domain in m
130 Notes
131 -----
132 To be used within sorted collections, instances of `_Subdomain` and its
133 subclasses need to be sortable. The order is defined with respect to the
134 left edge. Therefore, an equality test between two instances of
135 `_Subdomain` with the same `left_edge` attribute, and differing `length`
136 attributes, would hold. This unfortunately differs from the behaviour of
137 all other `attrs`-defined class, and can be surprising.
139 """
141 # TODO: total_ordering and SortedList are nice and all, but the surprise
142 # element of `==` not behaving as expected is not. Maybe consider a
143 # rewrite/perf comparison, doing without, and reinstatiating a list in
144 # order after breakup events, à la swisib.
146 left_edge: float
147 length: float
149 def __eq__(self, other: _Subdomain | Real) -> bool:
150 match other:
151 case _Subdomain():
152 return self.left_edge == other.left_edge
153 case Real(): 153 ↛ 155line 153 didn't jump to line 155, because the pattern on line 153 always matched
154 return self.left_edge == other
155 case _:
156 raise TypeError(
157 "Comparison not supported between instance of "
158 f"{type(self)} and {type(other)}"
159 )
161 def __lt__(self, other: _Subdomain | Real) -> bool:
162 match other:
163 case _Subdomain():
164 return self.left_edge < other.left_edge
165 case Real():
166 return self.left_edge < other
167 case _:
168 raise TypeError(
169 "Comparison not supported between instance of "
170 f"{type(self)} and {type(other)}"
171 )
173 @functools.cached_property
174 def right_edge(self):
175 return self.left_edge + self.length
178@attrs.define(frozen=True)
179class Ice:
180 """A container for ice mechanical properties.
182 Ice is modelled as an elastic thin plate, with prescribed density,
183 thickness, Poisson's ratio and Young's modulus. Its fracture under bending
184 is considered either through the lens of Griffith's fracture mechanics, or
185 through the framework of strain failure commonly used in the sea ice
186 modelling community. The fracture toughness is relevant to the former,
187 while the strain threshold is relevant to the latter. Ice is considered
188 translationally invariant in one horizontal direction, so that its
189 quadratic moment of area is given per unit length in that direction.
191 Parameters
192 ----------
193 density : float
194 Density in kg m^-3
195 frac_toughness : float
196 Fracture toughness in Pa m^-1/2
197 poissons_ratio : float
198 Poisson's ratio
199 strain_threshold : float
200 Critical flexural strain in m m^-1
201 thickness : float
202 Ice thickness in m
203 youngs_modulus : float
204 Scalar Young's modulus in Pa
206 Attributes
207 ----------
208 quad_moment : float
209 Quadratic moment of area in m^3
210 flex_rigidity : float
211 Flexural rigidity in
212 frac_energy_rate : float
213 Fracture energy release rate in J m^-2
215 """
217 density: float = 922.5
218 frac_toughness: float = 1e5
219 poissons_ratio: float = 0.3
220 strain_threshold: float = 3e-5
221 thickness: float = 1.0
222 youngs_modulus: float = 6e9
224 @functools.cached_property
225 def quad_moment(self) -> float:
226 return self.thickness**3 / (12 * (1 - self.poissons_ratio**2))
228 @functools.cached_property
229 def flex_rigidity(self) -> float:
230 return self.quad_moment * self.youngs_modulus
232 @functools.cached_property
233 def frac_energy_rate(self) -> float:
234 return (
235 (1 - self.poissons_ratio**2) * self.frac_toughness**2 / self.youngs_modulus
236 )
239@attrs.define(kw_only=True, frozen=True)
240class FloatingIce(Ice):
241 """An extension of `Ice` to represent properties due to buyoancy.
243 Parameters
244 ----------
245 draft : float
246 Immersed ice thickness at rest in m
247 dud : float
248 Height of the water column underneath the ice at rest in m
249 elastic_length_pow4 : float
250 Characteristic elastic length scale, raised to the 4th power, in m^4
252 Attributes
253 ----------
254 elastic_length : float
255 Characteristic elastic length scale in m
256 freeboard : float
257 Emerged ice thickness at rest in m
259 """
261 draft: float
262 dud: float
263 elastic_length_pow4: float
265 @classmethod
266 def from_ice_ocean(cls, ice: Ice, ocean: Ocean, gravity: float) -> FloatingIce:
267 """Build an instance by combining properties of existing objects.
269 Parameters
270 ----------
271 ice : Ice
272 ocean : Ocean
273 gravity : float
274 Strengh of the local gravitational field in m s^-2
276 Returns
277 -------
278 FloatingIce
280 """
281 draft = ice.density / ocean.density * ice.thickness
282 dud = ocean.depth - draft
283 # NOTE: as the 4th power of the elastic length scale arises naturally,
284 # we prefer using it to instantiate the class and computing the length
285 # scale when needed, over using the length scale for instantiation and
286 # recomputing the fourth power from it, as the latter approach can lead
287 # to substantial numerical imprecision.
288 el_lgth_pow4 = ice.flex_rigidity / (ocean.density * gravity)
289 return cls(
290 density=ice.density,
291 frac_toughness=ice.frac_toughness,
292 poissons_ratio=ice.poissons_ratio,
293 strain_threshold=ice.strain_threshold,
294 thickness=ice.thickness,
295 youngs_modulus=ice.youngs_modulus,
296 draft=draft,
297 dud=dud,
298 elastic_length_pow4=el_lgth_pow4,
299 )
301 @functools.cached_property
302 def elastic_length(self):
303 return self.elastic_length_pow4**0.25
305 @functools.cached_property
306 def freeboard(self):
307 return self.thickness - self.draft
309 @functools.cached_property
310 def _elastic_number(self) -> float:
311 """Reciprocal of the Characteristic elastic length scale.
313 Returns
314 -------
315 float
316 Elastic number in m^-1
318 """
319 return 1 / self.elastic_length
321 @functools.cached_property
322 def _red_elastic_number(self) -> float:
323 """Characteristic elastic number scaled by 1/sqrt(2).
325 Returns
326 -------
327 float
328 Reduced elastic number in m^-1
330 """
331 return 1 / (SQR2 * self.elastic_length)
334@attrs.define(frozen=True)
335class WavesUnderElasticPlate:
336 """A non-localised zone characterised by wave action.
338 The spatial behaviour of waves (wavelength) is linked to their temporal
339 behaviour (period) through a dispersion relation. In the case of waves
340 propagating underneath floating ice, considered as an elastic plate, this
341 dispersion relation depends on the properties of the ice as encapsulated by
342 the `FloatingIce` class.
344 Parameters
345 ----------
346 ice : FloatingIce
347 wavenumbers : 1d array_like of float
348 Propagating wavenumbers, in rad m^-1
350 """
352 ice: FloatingIce
353 wavenumbers: np.ndarray = attrs.field(repr=False)
355 @classmethod
356 def from_floating(
357 cls,
358 ice: FloatingIce,
359 spectrum: DiscreteSpectrum,
360 gravity: float,
361 ) -> typing.Self:
362 """Build an instance by combining properties of existing objects.
364 Parameters
365 ----------
366 ice : FloatingIce
367 spectrum : DiscreteSpectrum
368 gravity : float
369 Strengh of the local gravitational field in m s^-2
371 Returns
372 -------
373 WavesUnderElasticPlate
375 """
376 alphas = spectrum._ang_freqs_pow2 / gravity
377 deg1 = 1 - alphas * ice.draft
378 deg0 = -alphas * ice.elastic_length
379 scaled_ratio = ice.dud / ice.elastic_length
381 solver = dr.ElasticMassLoadingSolver(alphas, deg1, deg0, scaled_ratio)
382 # NOTE: `solver` returns non-dimensional wavenumbers
383 wavenumbers = solver.compute_wavenumbers() / ice.elastic_length
385 return cls(ice, wavenumbers)
387 @classmethod
388 def from_ocean(
389 cls,
390 ice: Ice,
391 ocean: Ocean,
392 spectrum: DiscreteSpectrum,
393 gravity: float,
394 ) -> typing.Self:
395 """Build an instance by combining properties of existing objects.
397 Parameters
398 ----------
399 ice : Ice
400 ocean : Ocean
401 spectrum : DiscreteSpectrum
402 gravity : float
403 Strengh of the local gravitational field in m s^-2
405 Returns
406 -------
407 WavesUnderElasticPlate
409 """
410 floating_ice = FloatingIce.from_ice_ocean(ice, ocean, gravity)
411 return cls.from_floating(floating_ice, spectrum, gravity)
414# TODO: check docstring
415@attrs.define(frozen=True)
416class WavesUnderIce:
417 """A non-localised zone characetrised by wave action under floating ice.
419 This class extends the behaviour of `WavesUnderElasticPlate` by adding an
420 `attenuations` attribute, that parametrises the observed exponential decay
421 of waves underneath floating ice.
423 Parameters
424 ----------
425 ice : FloatingIce
426 wavenumbers : 1d array_like of float
427 Propagating wavenumbers, in rad m^-1
428 attenuations : 1d array_like of float
429 Parametrised wave amplitude attenuation rate, in m^-1
431 """
433 ice: FloatingIce
434 wavenumbers: np.ndarray = attrs.field(repr=False)
435 attenuations: np.ndarray | Real = attrs.field(repr=False)
437 @classmethod
438 def without_attenuation(cls, waves_under_ep: WavesUnderElasticPlate) -> typing.Self:
439 """Build an instance by combining properties of existing objects.
441 Parameters
442 ----------
443 waves_under_ep : WavesUnderElasticPlate
444 An object instance
446 Returns
447 -------
448 WavesUnderIce
450 See Also
451 -----
452 lib.att.no_attenuation
454 """
455 return cls(
456 waves_under_ep.ice,
457 waves_under_ep.wavenumbers,
458 att.no_attenuation(),
459 )
461 @classmethod
462 def with_attenuation_01(cls, waves_under_ep: WavesUnderElasticPlate) -> typing.Self:
463 """Build an instance by combining properties of existing objects.
465 Parameters
466 ----------
467 waves_under_ep : WavesUnderElasticPlate
468 An object instance
470 Returns
471 -------
472 WavesUnderIce
474 See Also
475 -----
476 lib.att.parameterisation_01
478 """
479 return cls(
480 waves_under_ep.ice,
481 waves_under_ep.wavenumbers,
482 att.parameterisation_01(
483 waves_under_ep.ice.thickness, waves_under_ep.wavenumbers
484 ),
485 )
487 @classmethod
488 def with_generic_attenuation(
489 cls,
490 waves_under_ep: WavesUnderElasticPlate,
491 parameterisation: typing.Callable,
492 args: str | None = None,
493 **kwargs,
494 ) -> typing.Self:
495 """Instantiate a `WavesUnderFloe` with custom attenuation.
497 Parameters
498 ----------
499 waves_under_ep : WavesUnderElasticPlate
500 An object instance.
501 parameterisation : typing.Callable
502 Function defining attenuation.
503 Must return a type broadcastable to `waves_under_ep.wavenumbers`.
504 args : str | None
505 A string of attributes of `waves_under_ep`, separated by
506 whitespace, to be passed as parameters to `parameterisation`.
507 All parameters will be passed as a mapping between the stem of the
508 attribute and its value.
509 **kwargs : dict
510 Additional parameters to pass to `parameterisation`.
512 Returns
513 -------
514 WavesUnderFloe
516 Examples
517 --------
518 Assuming an existing `wue` instance of `WavesUnderElasticPlate`, the
519 three following objects are identical, setting the attenuation egal to
520 the ice density for all wave modes.
522 >>> WavesUnderIce.with_generic_attenuation_param(
523 wue,
524 lambda density: density,
525 "ice.density"
526 )
527 >>> WavesUnderIce.with_generic_attenuation_param(
528 wue,
529 lambda density: density,
530 {"density": wue.ice.density},
531 )
532 >>> WavesUnderIce(wue.ice, wue.wavenumbers, wue.ice.density)
534 """
535 if args is not None:
536 kwargs |= {
537 arg.split(".")[-1]: operator.attrgetter(arg)(waves_under_ep)
538 for arg in args.split()
539 }
540 return cls(
541 waves_under_ep.ice, waves_under_ep.wavenumbers, parameterisation(**kwargs)
542 )
544 @functools.cached_property
545 def _c_wavenumbers(self) -> np.ndarray:
546 """Complex wavenumbers.
548 Their real part correspond to the propagating wavenumber, while their
549 imaginary part correspond to the attenuation rate.
551 Returns
552 -------
553 1d np.ndarray of complex
554 The complex wavenumbers in m^-1
556 """
557 return self.wavenumbers + 1j * self.attenuations
560@attrs.define(frozen=True)
561class FreeSurfaceWaves:
562 """The wave state in the absence of ice.
564 The spatial behaviour of waves (wavelength) is linked to their temporal
565 behaviour (period) through a dispersion relation. In the case of free
566 surface waves, propagating underneath floating ice, this dispersion
567 relation depends on the properties of the ocean as encapsulated in the
568 `Ocean` class.
570 Parameters
571 ----------
572 ocean : Ocean
573 wavenumbers : array_like
574 Propagating wavenumbers in rad m^-1
576 Attributes
577 ----------
578 wavelengths : 1d np.ndarray of float
579 Propagating wavelengths in m
581 """
583 ocean: Ocean
584 wavenumbers: np.ndarray
586 @classmethod
587 def from_ocean(cls, ocean: Ocean, spectrum: DiscreteSpectrum, gravity: float):
588 """Build an instance by combining properties of existing objects."""
589 alphas = spectrum._ang_freqs_pow2 / gravity
590 solver = dr.FreeSurfaceSolver(alphas, ocean.depth)
591 wavenumbers = solver.compute_wavenumbers()
592 return cls(ocean, wavenumbers)
594 @functools.cached_property
595 def wavelengths(self) -> np.ndarray:
596 return PI_2 / self.wavenumbers
599# TODO: docstring inheritance
600@attrs.define(kw_only=True, eq=False)
601class Floe(_Subdomain):
602 """An ice floe localised in space.
604 Parameters
605 ----------
606 ice : Ice
607 The mechanical properties of the floe
609 """
611 ice: Ice
614@attrs.define(kw_only=True, eq=False)
615class WavesUnderFloe(_Subdomain):
616 """A localised zone characetrised by wave action under floating ice.
618 Parameters
619 ----------
620 wui : WavesUnderIce
621 edge_amplitudes : 1d np.ndarray of complex
622 The wave complex amplitude at the left edge of the floe in m
623 generation : int
624 The number of fractures that led to the existence of this floe
626 """
628 wui: WavesUnderIce
629 edge_amplitudes: np.ndarray
630 generation: int = 0
632 @functools.cached_property
633 def _adim(self) -> float:
634 """A non-dimentional number characetrising the floe.
636 Returns
637 -------
638 float
640 """
641 return self.length * self.wui.ice._red_elastic_number
643 # TODO: typing.Self?
644 def make_copy(self) -> WavesUnderFloe:
645 return WavesUnderFloe(
646 left_edge=self.left_edge,
647 length=self.length,
648 wui=self.wui,
649 edge_amplitudes=self.edge_amplitudes.copy(),
650 generation=self.generation,
651 )
653 @typing.overload
654 def shift_waves(self, phase_shifts: np.ndarray, inplace: typing.Literal[True]): ... 654 ↛ exitline 654 didn't return from function 'shift_waves'
656 # TODO: typing.Self
657 @typing.overload
658 def shift_waves( 658 ↛ exitline 658 didn't jump to the function exit
659 self, phase_shifts: np.ndarray, inplace: typing.Literal[False]
660 ) -> WavesUnderFloe: ...
662 # TODO: docstring
663 def shift_waves(self, phase_shifts: np.ndarray, inplace: bool = True):
664 shifted_amplitudes = self.edge_amplitudes * np.exp(-1j * phase_shifts)
665 if not inplace: 665 ↛ 666line 665 didn't jump to line 666, because the condition on line 665 was never true
666 return WavesUnderFloe(
667 left_edge=self.left_edge,
668 length=self.length,
669 wui=self.wui,
670 edge_amplitudes=shifted_amplitudes,
671 generation=self.generation,
672 )
673 # HACK: instantiating an new object is cheap, resorting a whole list of
674 # subdomains is not, so we mutate the amplitude/phase instead
675 object.__setattr__(self, "edge_amplitudes", shifted_amplitudes)
677 def displacement(
678 self,
679 x: np.ndarray,
680 growth_params=None,
681 an_sol: bool | None = None,
682 num_params=None,
683 ):
684 return ph.DisplacementHandler.from_wuf(self, growth_params).compute(
685 x, an_sol, num_params
686 )
688 def energy(self, growth_params=None, an_sol: bool | None = None, num_params=None):
689 return ph.EnergyHandler.from_wuf(self, growth_params).compute(
690 an_sol, num_params
691 )
694class DiscreteSpectrum:
695 def __init__(
696 self,
697 amplitudes,
698 frequencies,
699 phases=0,
700 ):
702 # np.ravel to force precisely 1D-arrays
703 # Promote the map to list so the iterator can be used several times
704 args = list(map(np.ravel, (amplitudes, frequencies, phases)))
705 (size,) = np.broadcast_shapes(*(arr.shape for arr in args))
707 # TODO: sort waves by frequencies or something
708 # TODO: sanity checks on nan, etc. that could be returned
709 # by the Spectrum objects
711 # If size is one, all the arguments are scalar and the "spectrum" is
712 # monochromatic. Otherwise, there is at least one argument with
713 # different components. The eventual other arguments with a single
714 # component are repeated for instantiating as many Wave objects as
715 # needed.
716 if size != 1:
717 for i, arr in enumerate(args):
718 if arr.size == 0: 718 ↛ 719line 718 didn't jump to line 719, because the condition on line 718 was never true
719 raise ValueError
720 if arr.size == 1:
721 args[i] = itertools.repeat(arr[0], size)
723 self.__waves = [
724 Wave.from_frequency(_a, _f, phase=_ph) for _a, _f, _ph in zip(*args)
725 ]
727 @property
728 def waves(self):
729 return self.__waves
731 @functools.cached_property
732 def _ang_freqs_pow2(self):
733 return self._ang_freqs**2
735 @functools.cached_property
736 def _amps(self):
737 return np.asarray([wave.amplitude for wave in self.waves])
739 @functools.cached_property
740 def _freqs(self):
741 return np.asarray([wave.frequency for wave in self.waves])
743 @functools.cached_property
744 def _ang_freqs(self):
745 return np.asarray([wave.angular_frequency for wave in self.waves])
747 @functools.cached_property
748 def _phases(self):
749 return np.asarray([wave.phase for wave in self.waves])
751 @functools.cached_property
752 def nf(self):
753 return len(self.waves)
756# TODO: docstrings
757@attrs.define
758class Domain:
759 """A spatial domain forced by waves.
761 This represents the state of a MIZ at a given time.
764 Attributes
765 ----------
766 gravity : float
767 spectrum : DiscreteSpectrum
768 fsw : FreeSurfaceWaves
769 attenuation: flexrac1d.lib.att.Attenuation
770 growth_params : list
771 subdomains : SortedList of WavesUnderFloe
772 cached_wuis :
773 cached_phases :
775 """
777 gravity: float
778 spectrum: DiscreteSpectrum
779 fsw: FreeSurfaceWaves
780 attenuation: att.Attenuation = attrs.field(repr=False)
781 growth_params: list[np.array, float] | None = None
782 subdomains: SortedList = attrs.field(repr=False, init=False, factory=SortedList)
783 cached_wuis: dict[Ice, WavesUnderIce] = attrs.field(
784 repr=False, init=False, factory=dict
785 )
786 cached_phases: dict[float, np.ndarray] = attrs.field(
787 repr=False, init=False, factory=dict
788 )
790 @classmethod
791 def from_discrete(
792 cls,
793 gravity,
794 spectrum,
795 ocean,
796 attenuation: att.Attenuation | None = None,
797 growth_params: tuple | None = None,
798 ):
799 fsw = FreeSurfaceWaves.from_ocean(ocean, spectrum, gravity)
800 if attenuation is None:
801 attenuation = att.AttenuationParameterisation(1)
802 return cls(gravity, spectrum, fsw, attenuation, growth_params)
804 @classmethod
805 def with_growth_means(
806 cls,
807 gravity: float,
808 spectrum: DiscreteSpectrum,
809 ocean: Ocean,
810 growth_means: np.ndarray | Sequence[Real] | Real,
811 attenuation: att.Attenuation | None = None,
812 ) -> typing.Self:
813 return cls.from_discrete(
814 gravity, spectrum, ocean, attenuation, (growth_means, None)
815 )
817 @classmethod
818 def with_growth_std(
819 cls,
820 gravity: float,
821 spectrum: DiscreteSpectrum,
822 ocean: Ocean,
823 growth_std: Real,
824 attenuation: att.Attenuation | None = None,
825 ) -> typing.Self:
826 return cls.from_discrete(gravity, spectrum, ocean, attenuation, (0, growth_std))
828 def __attrs_post_init__(self):
829 if self.growth_params is not None:
830 if len(self.growth_params) != 2:
831 raise ValueError
832 growth_means, growth_std = (
833 np.asarray(self.growth_params[0]),
834 self.growth_params[1],
835 )
836 # TODO: simplify all this. Ideally, do not test for anything or
837 # babysit the user. Why was upping growth_mean to a column
838 # necessary in case its of size 1?
839 if growth_means.size == 1:
840 # As `broadcast_to` returns a view,
841 # copying is necessary to obtain a mutable array. It is easier
842 # than dealing with 0-length and 1-length arrays seperately.
843 growth_means = np.broadcast_to(
844 growth_means, (self.spectrum.nf, 1)
845 ).copy()
846 else:
847 if growth_means.size != self.spectrum.nf: 847 ↛ 853line 847 didn't jump to line 853, because the condition on line 847 was never false
848 raise ValueError(
849 f"Means (size {growth_means.size}) could not be"
850 "broadcast with the shape of the spectrum"
851 f"({self.spectrum.nf} components)"
852 )
853 if growth_std is None:
854 growth_std = self.fsw.wavelengths[self.spectrum._amps.argmax()]
855 self.growth_params = [growth_means, growth_std]
857 def _compute_phase_shifts(self, delta_time: float):
858 if delta_time not in self.cached_phases:
859 self.cached_phases[delta_time] = delta_time * self.spectrum._ang_freqs
860 return self.cached_phases[delta_time]
862 def _compute_wui(self, ice: Ice):
863 if ice not in self.cached_wuis:
864 wup = WavesUnderElasticPlate.from_ocean(
865 ice, self.fsw.ocean, self.spectrum, self.gravity
866 )
867 if isinstance(self.attenuation, att.AttenuationParameterisation): 867 ↛ 873line 867 didn't jump to line 873, because the condition on line 867 was never false
868 if self.attenuation == att.AttenuationParameterisation.NO:
869 wui = WavesUnderIce.without_attenuation(wup)
870 elif self.attenuation == att.AttenuationParameterisation.PARAM_01: 870 ↛ 879line 870 didn't jump to line 879, because the condition on line 870 was never false
871 wui = WavesUnderIce.with_attenuation_01(wup)
872 else:
873 wui = WavesUnderIce.with_generic_attenuation(
874 wup,
875 self.attenuation.function,
876 self.attenuation.args,
877 **self.attenuation.kwargs,
878 )
879 self.cached_wuis[ice] = wui
880 return self.cached_wuis[ice]
882 def _shift_phases(self, phases: np.ndarray):
883 for i in range(len(self.floes)):
884 self.floes[i].phases -= phases
886 def _shift_growth_means(self, phases: np.ndarray):
887 # TODO: refine to take into account subdomain transitions
888 # and floes with variying properties
889 mask = self.growth_params[0] < self.subdomains[0].left_edge
890 if mask.any():
891 self.growth_params[0][mask] += (
892 phases[mask[:, 0]] / self.fsw.wavenumbers[mask[:, 0]]
893 )
894 if not mask.all():
895 self.growth_params[0][~mask] += (
896 phases[~mask[:, 0]] / self.subdomains[0].wui.wavenumbers[~mask[:, 0]]
897 )
899 def add_floes(self, floes: Floe | Sequence[Floe]):
900 subdomains = self._init_subdomains(floes)
901 self._add_c_floes(subdomains)
903 @staticmethod
904 def _promote_floe(floes: Floe | Sequence[Floe]) -> Sequence[Floe]:
905 match floes:
906 case Floe():
907 return (floes,)
908 case Sequence():
909 return floes
910 case _:
911 raise ValueError(
912 "`floes` should be a `Floe` object or a sequence of such objects"
913 )
915 def _check_overlap(self, floes: Sequence[Floe]):
916 l_edges, r_edges = map(
917 np.array, zip(*((floe.left_edge, floe.right_edge) for floe in floes))
918 )
919 if not (r_edges[:-1] <= l_edges[1:]).all():
920 raise ValueError("Floe overlap") # TODO: dedicated exception
922 def _init_phases(self, floes: Sequence[Floe]) -> np.ndarray:
923 phases = np.full((len(floes), self.spectrum.nf), np.nan)
924 phases[0] = self.spectrum._phases + floes[0].left_edge * self.fsw.wavenumbers
925 for i, floe in enumerate(floes[1:], 1):
926 wui = self._compute_wui(floe.ice)
927 prev = floes[i - 1]
928 phases[i:,] = (
929 phases[i - 1]
930 + floe.length * wui.wavenumbers
931 + (prev.right_edge - floe.left_edge) * self.fsw.wavenumbers
932 )
933 return phases % PI_2
935 def _init_amplitudes(self, floes: Sequence[Floe]) -> np.ndarray:
936 amplitudes = np.full((len(floes), self.spectrum.nf), np.nan)
937 amplitudes[0, :] = self.spectrum._amps
938 for i, floe in enumerate(floes[1:], 1):
939 amplitudes[i, :] = amplitudes[i - 1, :] * np.exp(
940 -self._compute_wui(floe.ice).attenuations * floe.length
941 )
942 return amplitudes
944 def _init_subdomains(self, floes: Sequence[Floe]) -> list[WavesUnderFloe]:
945 # TODO: look for already existing floes. In the present state, only
946 # valid for starting from scratch, not for adding floes to a domain
947 # that already has some.
948 floes = self.__class__._promote_floe(floes)
949 self._check_overlap(floes)
950 complex_amplitudes = self._init_amplitudes(floes) * np.exp(
951 1j * self._init_phases(floes)
952 )
954 return [
955 WavesUnderFloe(
956 left_edge=floe.left_edge,
957 length=floe.length,
958 wui=self._compute_wui(floe.ice),
959 edge_amplitudes=edge_amplitudes,
960 )
961 for floe, edge_amplitudes in zip(floes, complex_amplitudes)
962 ]
964 def iterate(self, delta_time: float):
965 phase_shifts = self._compute_phase_shifts(delta_time)
966 # TODO: can be optimised by iterating a first time to extract the
967 # edges, coerce them to a np.array, apply the product with
968 # complex_shifts, and then iterate a second time to build the objects.
969 # See Propagation_tests.ipynb/DNE06-26
970 for i in range(len(self.subdomains)):
971 self.subdomains[i].shift_waves(phase_shifts)
972 if self.growth_params is not None: 972 ↛ 974line 972 didn't jump to line 974, because the condition on line 972 was never true
973 # Phases are only modulo'd in the setter
974 self._shift_growth_means(phase_shifts)
976 def _pop_c_floe(self, wuf: WavesUnderFloe):
977 self.subdomains.remove(wuf)
979 def _add_c_floes(self, wuf: Sequence[WavesUnderFloe]):
980 # It is assume no overlap will occur, and phases have been properly
981 # set, as these method should only be called after a fracture event
982 self.subdomains.update(wuf)
984 def breakup(
985 self, fracture_handler: fh._FractureHandler, an_sol=None, num_params=None
986 ):
987 dct = {}
988 for i, wuf in enumerate(self.subdomains):
989 xf = fracture_handler.search(wuf, self.growth_params, an_sol, num_params)
990 if xf is not None:
991 old = wuf
992 new = fracture_handler.split(wuf, xf)
993 dct[i] = old, new
994 for old, new in dct.values():
995 self._pop_c_floe(old)
996 self._add_c_floes(new)
998 def plot(
999 self,
1000 resolution: float,
1001 left_bound: float,
1002 ax=None,
1003 an_sol=None,
1004 add_surface=True,
1005 base=0,
1006 kw_dis=None,
1007 kw_sur=None,
1008 ):
1009 plot_displacement(
1010 resolution, self, left_bound, ax, an_sol, add_surface, base, kw_dis, kw_sur
1011 )