Coverage for src/flexfrac1d/lib/displacement.py: 96%

1import numpy as np 

2from .constants import PI_D4, SQR2 

3from . import numerical 

4 

5 

6def _wavefield(x, c_amps, c_wavenumbers): 

7 return np.imag(c_amps @ _unit_wavefield(x, c_wavenumbers)) 

8 

9 

10def _unit_wavefield(x, c_wavenumbers): 

11 return np.exp((1j * c_wavenumbers[:, None]) * x) 

12 

13 

14def _dis_par_amps(red_num: float, wave_params: tuple[np.ndarray]): 

15 """Complex amplitudes of individual particular solutions""" 

16 c_amplitudes, c_wavenumbers = wave_params 

17 return c_amplitudes / (1 + 0.25 * (c_wavenumbers / red_num) ** 4) 

18 

19 

20def _dis_hom_mat(red_num: float, length: float): 

21 """Linear application to determine, from the BCs, the coefficients of 

22 the four independent solutions to the homo ODE""" 

23 adim = red_num * length 

24 adim2 = 2 * adim 

25 denom = ( 

26 -2 

27 * red_num**2 

28 * (np.expm1(-adim2) ** 2 + 2 * np.exp(-adim2) * (np.cos(adim2) - 1)) 

29 ) 

30 mat = np.array( 

31 [ 

32 [ 

33 SQR2 * np.exp(-adim) * np.sin(adim2 + PI_D4) - np.exp(-3 * adim), 

34 -SQR2 * np.cos(adim + PI_D4) 

35 - 3 * np.exp(-adim2) * np.sin(adim) 

36 + np.exp(-adim2) * np.cos(adim), 

37 np.exp(-adim) * np.sin(adim2) - np.exp(-adim) + np.exp(-3 * adim), 

38 np.cos(adim) 

39 - 2 * np.exp(-adim2) * np.sin(adim) 

40 - np.exp(-adim2) * np.cos(adim), 

41 ], 

42 [ 

43 -SQR2 * np.exp(-adim) * np.cos(adim2 + PI_D4) 

44 + 2 * np.exp(-adim) 

45 - np.exp(-3 * adim), 

46 -SQR2 * np.sin(adim + PI_D4) 

47 + SQR2 * np.exp(-adim2) * np.cos(adim + PI_D4), 

48 -np.exp(-adim) * np.cos(adim2) + np.exp(-adim), 

49 np.sin(adim) - np.exp(-adim2) * np.sin(adim), 

50 ], 

51 [ 

52 -1 + SQR2 * np.exp(-adim2) * np.cos(adim2 + PI_D4), 

53 (3 * np.sin(adim) + np.cos(adim)) * np.exp(-adim) 

54 - SQR2 * np.exp(-3 * adim) * np.sin(adim + PI_D4), 

55 (np.sin(adim2) + 1) * np.exp(-adim2) - 1, 

56 (-2 * np.sin(adim) + np.cos(adim)) * np.exp(-adim) 

57 - np.exp(-3 * adim) * np.cos(adim), 

58 ], 

59 [ 

60 (SQR2 * np.sin(adim2 + PI_D4) - 2) * np.exp(-adim2) + 1, 

61 -SQR2 * np.exp(-adim) * np.sin(adim + PI_D4) 

62 + SQR2 * np.exp(-3 * adim) * np.cos(adim + PI_D4), 

63 (-np.cos(adim2) + 1) * np.exp(-adim2), 

64 np.exp(-adim) * np.sin(adim) - np.exp(-3 * adim) * np.sin(adim), 

65 ], 

66 ] 

67 ) 

68 mat[:, 2:] /= red_num 

69 

70 return mat / denom 

71 

72 

73def _dis_hom_rhs(floe_params: tuple[float], wave_params: tuple[np.ndarray]): 

74 """Vector onto which apply the matrix, to extract the coefficients""" 

75 red_num, length = floe_params 

76 _, c_wavenumbers = wave_params 

77 exp_arg = 1j * c_wavenumbers * length 

78 

79 r1 = c_wavenumbers**2 * _dis_par_amps(red_num, wave_params) 

80 r2 = np.imag(r1 @ np.exp(exp_arg)) 

81 r3 = 1j * c_wavenumbers * r1 

82 r4 = np.imag(r3 @ np.exp(exp_arg)) 

83 r1 = np.imag(r1).sum() 

84 r3 = np.imag(r3).sum() 

85 

86 return np.array((r1, r2, r3, r4)) 

87 

88 

89def _dis_hom_coefs( 

90 floe_params: tuple[float], wave_params: tuple[np.ndarray] 

91) -> np.ndarray: 

92 """Coefficients of the four orthogonal homogeneous solutions""" 

93 return _dis_hom_mat(*floe_params) @ _dis_hom_rhs(floe_params, wave_params) 

94 

95 

96def _dis_hom(x: np.ndarray, floe_params: tuple[float], wave_params: tuple[np.ndarray]): 

97 """Homogeneous solution to the displacement ODE""" 

98 red_num, length = floe_params 

99 arr = red_num * x 

100 cosx, sinx = np.cos(arr), np.sin(arr) 

101 expx = np.exp(-red_num * (length - x)) 

102 exmx = np.exp(-arr) 

103 return np.vstack( 

104 ([expx * cosx, expx * sinx, exmx * cosx, exmx * sinx]) 

105 ).T @ _dis_hom_coefs(floe_params, wave_params) 

106 

107 

108def _dis_par(x: np.ndarray, red_num: float, wave_params: tuple[np.ndarray]): 

109 """Sum of the particular solutions to the displacement ODE""" 

110 return _wavefield(x, _dis_par_amps(red_num, wave_params), wave_params[1]) 

111 

112 

113def displacement( 

114 x, 

115 floe_params: tuple[float], 

116 wave_params: tuple[np.ndarray], 

117 growth_params: tuple | None = None, 

118 an_sol: bool | None = None, 

119 num_params: dict | None = None, 

120): 

121 if numerical._use_an_sol(an_sol, floe_params[1], growth_params): 121 ↛ 125line 121 didn't jump to line 125, because the condition on line 121 was never false

122 return _dis_hom(x, floe_params, wave_params) + _dis_par( 

123 x, floe_params[0], wave_params 

124 ) 

125 return numerical.displacement( 

126 x, floe_params, wave_params, growth_params, num_params 

127 )