Coverage for src/flexfrac1d/lib/physics.py: 97%
228 statements
« prev ^ index » next coverage.py v7.4.1, created at 2024-08-30 14:00 +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(wuf: model.WavesUnderFloe, growth_params):
18 floe_params = wuf.wui.ice._red_elastic_number, wuf.length
19 wave_params = wuf.edge_amplitudes, wuf.wui._c_wavenumbers
20 if growth_params is not None: 20 ↛ 21line 20 didn't jump to line 21, because the condition on line 20 was never true
21 growth_params = growth_params[0] - wuf.left_edge, growth_params[1]
22 return floe_params, wave_params, growth_params
25def _dis_par_amps(red_num: float, wave_params: tuple[np.ndarray]):
26 """Complex amplitudes of individual particular solutions"""
27 c_amplitudes, c_wavenumbers = wave_params
28 return c_amplitudes / (1 + 0.25 * (c_wavenumbers / red_num) ** 4)
31def _dis_hom_coefs(
32 floe_params: tuple[float], wave_params: tuple[np.ndarray]
33) -> np.ndarray:
34 """Coefficients of the four orthogonal homogeneous solutions"""
35 return _dis_hom_mat(*floe_params) @ _dis_hom_rhs(floe_params, wave_params)
38def _dis_hom_mat(red_num: float, length: float):
39 """Linear application to determine, from the BCs, the coefficients of
40 the four independent solutions to the homo ODE"""
41 adim = red_num * length
42 adim2 = 2 * adim
43 denom = (
44 -2
45 * red_num**2
46 * (np.expm1(-adim2) ** 2 + 2 * np.exp(-adim2) * (np.cos(adim2) - 1))
47 )
48 mat = np.array(
49 [
50 [
51 SQR2 * np.exp(-adim) * np.sin(adim2 + PI_D4) - np.exp(-3 * adim),
52 -SQR2 * np.cos(adim + PI_D4)
53 - 3 * np.exp(-adim2) * np.sin(adim)
54 + np.exp(-adim2) * np.cos(adim),
55 np.exp(-adim) * np.sin(adim2) - np.exp(-adim) + np.exp(-3 * adim),
56 np.cos(adim)
57 - 2 * np.exp(-adim2) * np.sin(adim)
58 - np.exp(-adim2) * np.cos(adim),
59 ],
60 [
61 -SQR2 * np.exp(-adim) * np.cos(adim2 + PI_D4)
62 + 2 * np.exp(-adim)
63 - np.exp(-3 * adim),
64 -SQR2 * np.sin(adim + PI_D4)
65 + SQR2 * np.exp(-adim2) * np.cos(adim + PI_D4),
66 -np.exp(-adim) * np.cos(adim2) + np.exp(-adim),
67 np.sin(adim) - np.exp(-adim2) * np.sin(adim),
68 ],
69 [
70 -1 + SQR2 * np.exp(-adim2) * np.cos(adim2 + PI_D4),
71 (3 * np.sin(adim) + np.cos(adim)) * np.exp(-adim)
72 - SQR2 * np.exp(-3 * adim) * np.sin(adim + PI_D4),
73 (np.sin(adim2) + 1) * np.exp(-adim2) - 1,
74 (-2 * np.sin(adim) + np.cos(adim)) * np.exp(-adim)
75 - np.exp(-3 * adim) * np.cos(adim),
76 ],
77 [
78 (SQR2 * np.sin(adim2 + PI_D4) - 2) * np.exp(-adim2) + 1,
79 -SQR2 * np.exp(-adim) * np.sin(adim + PI_D4)
80 + SQR2 * np.exp(-3 * adim) * np.cos(adim + PI_D4),
81 (-np.cos(adim2) + 1) * np.exp(-adim2),
82 np.exp(-adim) * np.sin(adim) - np.exp(-3 * adim) * np.sin(adim),
83 ],
84 ]
85 )
86 mat[:, 2:] /= red_num
88 return mat / denom
91def _dis_hom_rhs(floe_params: tuple[float], wave_params: tuple[np.ndarray]):
92 """Vector onto which apply the matrix, to extract the coefficients"""
93 red_num, length = floe_params
94 _, c_wavenumbers = wave_params
95 exp_arg = 1j * c_wavenumbers * length
97 r1 = c_wavenumbers**2 * _dis_par_amps(red_num, wave_params)
98 r2 = np.imag(r1 @ np.exp(exp_arg))
99 r3 = 1j * c_wavenumbers * r1
100 r4 = np.imag(r3 @ np.exp(exp_arg))
101 r1 = np.imag(r1).sum()
102 r3 = np.imag(r3).sum()
104 return np.array((r1, r2, r3, r4))
107# Can be used to decorate methods to make them return a scalar instead of
108# a 1d-array of length 1, if the argument is a scalar OR a 0d array (that is, a
109# scalar). If the argument is a tuple or a list or a 1d array of length 1,
110# a 1d-array of length 1 is returned.
111# This emulate the behaviour of standard numpy ufuncs.
112def _demote_to_scalar(f):
113 @functools.wraps(f)
114 def wrapper(*args):
115 instance, x = args
116 x = np.asarray(x)
117 res = f(instance, x)
118 if len(res) == 1 and x.ndim == 0:
119 return res.item()
120 return res
122 return wrapper
125@attrs.define
126class FluidSurfaceHandler:
127 wave_params: tuple[np.ndarray]
128 growth_params: list[np.ndarray, float] | None = None
130 @classmethod
131 def from_wuf(cls, wuf: model.WavesUnderFloe, growth_params=None):
132 return cls(*(_package_wuf(wuf, growth_params)[1:]))
134 @classmethod
135 def from_domain(cls, domain: model.Domain, growth_params: list | None = None):
136 return cls(
137 (
138 domain.spectrum._amps * np.exp(1j * domain.spectrum._phases),
139 domain.fsw.wavenumbers,
140 ),
141 growth_params,
142 )
144 @_demote_to_scalar
145 def compute(self, x, /):
146 return numerical.free_surface(x, self.wave_params, self.growth_params)
149@attrs.define
150class DisplacementHandler:
151 floe_params: tuple[float]
152 wave_params: tuple[np.ndarray]
153 growth_params: list[np.ndarray, float] | None = None
155 @classmethod
156 def from_wuf(cls, wuf: model.WavesUnderFloe, growth_params=None):
157 return cls(*_package_wuf(wuf, growth_params))
159 def _dis_hom(self, x: np.ndarray):
160 """Homogeneous solution to the displacement ODE"""
161 red_num, length = self.floe_params
162 arr = red_num * x
163 cosx, sinx = np.cos(arr), np.sin(arr)
164 expx = np.exp(-red_num * (length - x))
165 exmx = np.exp(-arr)
166 return np.vstack(
167 ([expx * cosx, expx * sinx, exmx * cosx, exmx * sinx])
168 ).T @ _dis_hom_coefs(self.floe_params, self.wave_params)
170 def _dis_par(self, x: np.ndarray):
171 """Sum of the particular solutions to the displacement ODE"""
172 return _wavefield(
173 x, _dis_par_amps(self.floe_params[0], self.wave_params), self.wave_params[1]
174 )
176 @_demote_to_scalar
177 def _dis(self, x, /):
178 return self._dis_hom(x) + self._dis_par(x)
180 def compute(
181 self,
182 x,
183 an_sol: bool | None = None,
184 num_params: dict | None = None,
185 ):
186 if numerical._use_an_sol(an_sol, self.floe_params[1], self.growth_params):
187 return self._dis(np.asarray(x))
188 return numerical.displacement(
189 x, self.floe_params, self.wave_params, self.growth_params, num_params
190 )
193@attrs.define
194class CurvatureHandler:
195 floe_params: tuple[float]
196 wave_params: tuple[np.ndarray]
197 growth_params: list[np.ndarray, float] | None = None
199 @classmethod
200 def from_wuf(cls, wuf: model.WavesUnderFloe, growth_params=None):
201 return cls(*_package_wuf(wuf, growth_params))
203 def _cur_wavefield(self, x):
204 """Second derivative of the interface"""
205 red_num = self.floe_params[0]
206 _, c_wavenumbers = self.wave_params
208 return -np.imag(
209 (_dis_par_amps(red_num, self.wave_params) * c_wavenumbers**2)
210 @ np.exp(1j * c_wavenumbers[:, None] * x)
211 )
213 def _cur_hom(self, x):
214 """Second derivative of the homogeneous part of the displacement"""
215 red_num, length = self.floe_params
216 arr = red_num * x
217 cosx, sinx = np.cos(arr), np.sin(arr)
218 expx = np.exp(-red_num * (length - x))
219 exmx = np.exp(-arr)
220 return (
221 2
222 * red_num**2
223 * np.vstack(([-expx * sinx, expx * cosx, exmx * sinx, -exmx * cosx])).T
224 @ _dis_hom_coefs(self.floe_params, self.wave_params)
225 )
227 def _cur_par(self, x):
228 """Second derivative of the particular part of the displacement"""
229 return self._cur_wavefield(x)
231 @_demote_to_scalar
232 def _cur(self, x, /):
233 return self._cur_hom(x) + self._cur_par(x)
235 def compute(
236 self,
237 x,
238 an_sol: bool | None = None,
239 num_params: dict | None = None,
240 ):
241 """Curvature of the floe, i.e. second derivative of the vertical displacement"""
242 if numerical._use_an_sol(an_sol, self.floe_params[1], self.growth_params):
243 return self._cur(x)
244 return numerical.curvature(
245 x, self.floe_params, self.wave_params, self.growth_params, num_params
246 )
249@attrs.define
250class StrainHandler:
251 curv_handler: CurvatureHandler
252 thickness: float
254 @classmethod
255 def from_wuf(cls, wuf: model.WavesUnderFloe, growth_params=None):
256 return cls(CurvatureHandler.from_wuf(wuf, growth_params), wuf.wui.ice.thickness)
258 def compute(self, x, an_sol: bool | None = None, num_params: dict | None = None):
259 return -self.thickness / 2 * self.curv_handler.compute(x, an_sol, num_params)
262@attrs.define
263class EnergyHandler:
264 floe_params: tuple[float]
265 wave_params: tuple[np.ndarray]
266 growth_params: list[np.ndarray, float] | None = None
267 factor: float = 1
269 @classmethod
270 def from_wuf(cls, wuf: model.WavesUnderFloe, growth_params=None):
271 factor = wuf.wui.ice.flex_rigidity / (2 * wuf.wui.ice.thickness)
272 return cls(*_package_wuf(wuf, growth_params), factor)
274 def _egy_hom(self):
275 """Energy from the homogen term of the displacement ODE"""
276 red_num, length = self.floe_params
277 adim = red_num * length
278 adim2 = 2 * adim
279 c_1, c_2, c_3, c_4 = _dis_hom_coefs(self.floe_params, self.wave_params)
281 return red_num**3 * (
282 +(c_1**2) * (-SQR2 * np.sin(adim2 + PI_D4) + 2 - np.exp(-adim2)) / 2
283 + c_2**2 * (SQR2 * np.sin(adim2 + PI_D4) + 2 - 3 * np.exp(-adim2)) / 2
284 + c_3**2 * ((SQR2 * np.cos(adim2 + PI_D4) - 2) * np.exp(-adim2) + 1) / 2
285 + c_4**2 * (3 - (SQR2 * np.cos(adim2 + PI_D4) + 2) * np.exp(-adim2)) / 2
286 + c_1 * c_2 * (SQR2 * np.cos(adim2 + PI_D4) - np.exp(-adim2))
287 - 2
288 * red_num
289 * c_1
290 * c_3
291 * (2 * length - np.sin(adim2) / red_num)
292 * np.exp(-adim)
293 + 4 * c_1 * c_4 * np.exp(-adim) * np.sin(adim) ** 2
294 + 4 * c_2 * c_3 * np.exp(-adim) * np.sin(adim) ** 2
295 - 2
296 * red_num
297 * c_2
298 * c_4
299 * (2 * length + np.sin(adim2) / red_num)
300 * np.exp(-adim)
301 + c_3 * c_4 * (-1 + SQR2 * np.exp(-adim2) * np.sin(adim2 + PI_D4))
302 )
304 def _egy_par_vals(self):
305 red_num = self.floe_params[0]
306 comp_amps = _dis_par_amps(red_num, self.wave_params)
307 _, c_wavenumbers = self.wave_params
309 return comp_amps * (1j * c_wavenumbers) ** 2
311 def _egy_par_pow2(self):
312 """Energy contribution from individual forcings"""
313 red_num, length = self.floe_params
314 _, c_wavenumbers = self.wave_params
315 wavenumbers = np.real(c_wavenumbers)
316 attenuations = np.imag(c_wavenumbers)
318 comp_curvs = self._egy_par_vals()
319 wn_moduli, curv_moduli = map(np.abs, (c_wavenumbers, comp_curvs))
320 wn_phases, curv_phases = map(np.angle, (c_wavenumbers, comp_curvs))
322 # Because of a division by `attenuations`, this term must me handled
323 # with some extra care
324 attenuation_mask = attenuations != 0
325 non_oscillatory_term = np.full(attenuations.shape, np.nan)
326 non_oscillatory_term[attenuation_mask] = (
327 -np.expm1(-2 * attenuations[attenuation_mask] * length)
328 / attenuations[attenuation_mask]
329 )
330 # NOTE: corresponds to the limit of the previous term, when attenuation -> 0
331 non_oscillatory_term[~attenuation_mask] = 2 * length
333 return (
334 curv_moduli**2
335 @ (
336 non_oscillatory_term
337 + (
338 np.sin(2 * curv_phases - wn_phases)
339 - np.sin(2 * (wavenumbers * length + curv_phases) - wn_phases)
340 * np.exp(-2 * attenuations * length)
341 )
342 / wn_moduli
343 )
344 ) / 4
346 def _egy_par_m(self):
347 """Energy contribution from forcing interactions"""
348 red_num, length = self.floe_params
349 _, c_wavenumbers = self.wave_params
350 wavenumbers = np.real(c_wavenumbers)
351 attenuations = np.imag(c_wavenumbers)
352 comp_curvs = self._egy_par_vals()
354 # Binomial coefficients, much quicker than itertools
355 idx1, idx2 = np.triu_indices(c_wavenumbers.size, 1)
357 mean_attenuations = attenuations[idx1] + attenuations[idx2]
358 comp_wns = (
359 wavenumbers[idx1] - wavenumbers[idx2],
360 wavenumbers[idx1] + wavenumbers[idx2],
361 )
362 comp_wns += 1j * mean_attenuations
364 curv_moduli = np.abs(comp_curvs)
365 _curv_phases = np.angle(comp_curvs)
366 curv_phases = (
367 _curv_phases[idx1] - _curv_phases[idx2],
368 _curv_phases[idx1] + _curv_phases[idx2],
369 )
371 def _f(comp_wns, curv_phases):
372 wn_moduli = np.abs(comp_wns)
373 wn_phases = np.angle(comp_wns)
374 return (
375 np.sin(curv_phases - wn_phases)
376 - (
377 np.exp(-np.imag(comp_wns) * length)
378 * np.sin(np.real(comp_wns) * length + curv_phases - wn_phases)
379 )
380 ) / wn_moduli
382 return (
383 curv_moduli[idx1]
384 * curv_moduli[idx2]
385 @ (_f(comp_wns[1], curv_phases[1]) - _f(comp_wns[0], curv_phases[0]))
386 )
388 def _egy_par(self):
389 """Energy from the forcing term of the displacement ODE"""
390 return self._egy_par_pow2() + self._egy_par_m()
392 def _egy_m(self):
393 red_num, length = self.floe_params
394 _, c_wavenumbers = self.wave_params
395 adim = red_num * length
396 wavenumbers = np.real(c_wavenumbers)
397 attenuations = np.imag(c_wavenumbers)
399 wn_abs, wn_phases = (_f(c_wavenumbers) for _f in (np.abs, np.angle))
400 comp_curvs = self._egy_par_vals()
401 curv_moduli, curv_phases = np.abs(comp_curvs), np.angle(comp_curvs)
402 num_add = red_num + wavenumbers
403 num_sub = red_num - wavenumbers
404 att_add = red_num + attenuations
405 att_sub = red_num - attenuations
406 q_pp = att_add**2 + num_add**2
407 q_mp = att_sub**2 + num_add**2
408 q_pm = att_add**2 + num_sub**2
409 q_mm = att_sub**2 + num_sub**2
411 phase_deltas = wn_phases - curv_phases
412 sin_phase_deltas = np.sin(phase_deltas)
413 cos_phase_deltas = np.cos(phase_deltas)
414 phase_quads = curv_phases + PI_D4
415 sin_phase_quads = np.sin(phase_quads)
416 cos_phase_quads = np.cos(phase_quads)
417 energ_1_K = np.array(
418 (
419 sin_phase_deltas * (1 / q_mp - 1 / q_mm),
420 cos_phase_deltas * (1 / q_mp + 1 / q_mm),
421 -sin_phase_deltas * (1 / q_pp - 1 / q_pm),
422 -cos_phase_deltas * (1 / q_pp + 1 / q_pm),
423 )
424 )
425 energ_1_b = np.array(
426 (
427 -(sin_phase_quads / q_mp - cos_phase_quads / q_mm),
428 (cos_phase_quads / q_mp - sin_phase_quads / q_mm),
429 -(cos_phase_quads / q_pp - sin_phase_quads / q_pm),
430 -(sin_phase_quads / q_pp - cos_phase_quads / q_pm),
431 )
432 )
433 energ_1 = wn_abs * energ_1_K + SQR2 * red_num * energ_1_b
434 energ_1[:2, :] *= np.exp(-adim)
436 arg_add = num_add * length
437 arg_add_phase = arg_add - phase_deltas
438 sin_arg_ap = np.sin(arg_add_phase)
439 cos_arg_ap = np.cos(arg_add_phase)
440 arg_sub = num_sub * length
441 arg_sub_phase = arg_sub + phase_deltas
442 sin_arg_sp = np.sin(arg_sub_phase)
443 cos_arg_sp = np.cos(arg_sub_phase)
444 arg_add_quad = arg_add + phase_quads
445 sin_arg_aq = np.sin(arg_add_quad)
446 cos_arg_aq = np.cos(arg_add_quad)
447 arg_sub_quad = arg_sub - curv_phases + PI_D4
448 sin_arg_sq = np.sin(arg_sub_quad)
449 cos_arg_sq = np.cos(arg_sub_quad)
450 energ_exp_K = np.array(
451 (
452 sin_arg_ap / q_mp + sin_arg_sp / q_mm,
453 -(cos_arg_ap / q_mp + cos_arg_sp / q_mm),
454 -(sin_arg_ap / q_pp + sin_arg_sp / q_pm),
455 cos_arg_ap / q_pp + cos_arg_sp / q_pm,
456 )
457 )
458 energ_exp_b = np.array(
459 (
460 sin_arg_aq / q_mp - sin_arg_sq / q_mm,
461 -cos_arg_aq / q_mp + cos_arg_sq / q_mm,
462 cos_arg_aq / q_pp - cos_arg_sq / q_pm,
463 sin_arg_aq / q_pp - sin_arg_sq / q_pm,
464 )
465 )
466 energ_exp = wn_abs * energ_exp_K + SQR2 * red_num * energ_exp_b
467 energ_exp[2:, :] *= np.exp(-adim)
468 energ_exp *= np.exp(-attenuations * length)
470 energ = (
471 (energ_1 + energ_exp)
472 @ curv_moduli
473 @ _dis_hom_coefs(self.floe_params, self.wave_params)
474 )
475 return red_num**2 * energ
477 def compute(
478 self,
479 an_sol: bool | None = None,
480 num_params: dict | None = None,
481 ):
482 if numerical._use_an_sol(an_sol, self.floe_params[1], self.growth_params): 482 ↛ 485line 482 didn't jump to line 485, because the condition on line 482 was never false
483 unit_energy = self._egy_hom() + 2 * self._egy_m() + self._egy_par()
484 else:
485 unit_energy = numerical.energy(
486 self.floe_params, self.wave_params, self.growth_params, num_params
487 )[0]
489 return self.factor * unit_energy