To receive notifications about scheduled maintenance, please subscribe to the mailing-list gitlab-operations@sympa.ethz.ch. You can subscribe to the mailing-list at https://sympa.ethz.ch

Commit 44c6df98 authored by Cyrill Mast's avatar Cyrill Mast
Browse files

Bugfix: issue #29

parent e6aeb072
Pipeline #88252 failed with stages
in 2 minutes and 18 seconds
......@@ -77,20 +77,24 @@ class ARF(BaseYosioka, BaseARF, BaseSphereFrequencyComposite):
:rtype: float
"""
if self.particle_type == "analytical":
return self._compute_analytical_arf()
if self.particle_type == "numerical":
return self._compute_numerical_arf()
if self.particle_type == "special":
return self._compute_special_arf()
out = self._compute_analytical_arf()
elif self.particle_type == "numerical":
out = self._compute_numerical_arf()
elif self.particle_type == "special":
out = self._compute_special_arf()
elif self.particle_type == "bubble":
return self._compute_bubble_arf()
out = self._compute_bubble_arf()
elif self.particle_type == "sphere":
return self._compute_sphere_arf()
out = self._compute_sphere_arf()
else:
print(self.particle_type)
raise NotImplementedError('Solution is no implemented (yet).')
surface = 4 * pi * self.R_0 ** 2
return out * surface
def _compute_analytical_arf(self) -> float:
"""Computes general analytical solution for the ARF according to
eq. (44)
......@@ -98,43 +102,41 @@ class ARF(BaseYosioka, BaseARF, BaseSphereFrequencyComposite):
"""
if self.n_max is None:
self.n_max = 10
self.n = 0
n = 0
out = [0 for _ in range(5)]
while self.n < self.n_max:
K_n_prev = self.K_n_0
M_n_prev = self.M_n_0
self.n += 1
while n < self.n_max:
# (45)
s1 = 0.5 * (K_n_prev * conj(self.K_n_0)).real
s1 *= 2 * self.n
s1 *= (-2 * pi * self.R_0 ** 2 * self.rho_f)
s1 = (self.K_n(n) * conj(self.K_n(n + 1))).real
s1 *= self.n + 1
out[0] += s1
# (46)
s2 = 0.5 * (M_n_prev * conj(self.M_n_0)).real
s2 *= 2 * (self.n - 1) * self.n * (self.n + 1)
s2 *= (2 * pi * self.xlambda ** 2 * self.rho_f)
s2 = (self.M_n(n) * conj(self.M_n(n + 1))).real
s2 *= n * (n + 1) * (n + 2)
out[1] += s2
# (47)
s3 = 0.5 * (K_n_prev * conj(self.M_n_0)).real
s4 = 0.5 * (M_n_prev * conj(self.K_n_0)).real
s3 *= 2 * self.n * (self.n + 1)
s3 *= (-2 * pi * self.R_0 * self.xlambda * self.rho_f)
s4 *= 2 * (self.n - 1) * self.n
s4 *= 2 * pi * self.R_0 * self.xlambda * self.rho_f
s3 = (self.K_n(n) * conj(self.M_n(n + 1))).real
s3 *= (n + 1) * (n + 2)
s4 = (self.M_n(n) * conj(self.K_n(n + 1))).real
s4 *= n * (n + 1)
out[2] += s3
out[3] += s4
# (48)
s5 = 0.5 * (M_n_prev * conj(self.M_n_0)).real
s5 *= 2 * self.n
s5 *= (-2 * pi * (self.k_f ** 2) *
(self.R_0 ** 2) * (self.xlambda ** 2) * self.rho_f)
s5 = (self.M_n(n) * conj(self.M_n(n + 1))).real
s5 *= n + 1
out[4] += s5
# print(out)
n += 1
out[0] *= -2 * pi * self.R_0 ** 2 * self.rho_f
out[1] *= 2 * pi * self.xlambda ** 2 * self.rho_f
out[2] *= -2 * pi * self.R_0 * self.xlambda * self.rho_f
out[3] *= 2 * pi * self.R_0 * self.xlambda * self.rho_f
out[4] *= (-2 * pi * (self.k_f ** 2) *
(self.R_0 ** 2) * (self.xlambda ** 2) * self.rho_f)
# (44)
return sum(out)
......@@ -285,8 +287,6 @@ class ARF(BaseYosioka, BaseARF, BaseSphereFrequencyComposite):
text = 'Criteria for small spheres might not be fulfilled.'
text += 'Consider using another model approach'
log.info(text)
print(text)
print(i)
break
if self.wave_type == 'standing':
return self._standing_wave_solution_sphere()
......@@ -321,8 +321,6 @@ class ARF(BaseYosioka, BaseARF, BaseSphereFrequencyComposite):
text = 'Criteria for small bubbles might not be fulfilled.'
text += 'Consider using another model approach'
log.info(text)
print(text)
print(i)
break
if self.wave_type == 'standing':
......
......@@ -2,7 +2,7 @@ from __future__ import annotations
from typing import Union
from gorkov import (
log, pi, exp,
log, pi, exp, full_range,
BaseSphereFrequencyComposite,
Sphere, Frequency,
BackgroundField,
......@@ -10,6 +10,7 @@ from gorkov import (
ActiveVariable, PassiveVariable,
StringFormatter as SF,
SpecialFunctions as sp
)
......@@ -57,7 +58,6 @@ class BaseYosioka(BaseSphereFrequencyComposite):
# Independent variables
self.n_max = n_max
self.n = 0
self.h = 0
self._particle_type = PassiveVariable(
particle_type,
......@@ -67,7 +67,6 @@ class BaseYosioka(BaseSphereFrequencyComposite):
self._n = PassiveVariable(self.n, "Index of iteration")
self._r = PassiveVariable(position, "Position")
self._omega = PassiveVariable(2 * pi * frequency, "Angular Frequency")
self._h = PassiveVariable(0, "height in standing field")
self._t = PassiveVariable(0, "time")
self._R_0 = PassiveVariable(radius, "Particle-Radius")
# Dependent variables
......@@ -77,20 +76,20 @@ class BaseYosioka(BaseSphereFrequencyComposite):
"ratio of speed of sounds")
self._F = ActiveVariable(self._compute_F,
"density-compressibility factor")
self._A_n = ActiveVariable(self._compute_A_n,
self._A_n = ActiveVariable(self._reset_A_n,
"Coefficient A_n for traveling field")
self._B_n = ActiveVariable(self._compute_B_n,
self._B_n = ActiveVariable(self._reset_B_n,
"Coefficient B_n for traveling field")
self._C_n = ActiveVariable(self._compute_C_n,
self._C_n = ActiveVariable(self._reset_C_n,
"Coefficient C_n for standing field")
self._D_n = ActiveVariable(self._compute_D_n,
self._D_n = ActiveVariable(self._reset_D_n,
"Coefficient C_n for standing field")
self._delta_n = ActiveVariable(self._compute_delta_n,
self._delta_n = ActiveVariable(self._reset_delta_n,
"Coefficient delta_n for calculations")
self._M_n_0 = ActiveVariable(self._compute_M_n_0,
"Coefficient M_n_0")
self._K_n_0 = ActiveVariable(self._compute_K_n_0,
"Coefficient K_n_0")
self._M_n = ActiveVariable(self._reset_M_n,
"Coefficient M_n_0")
self._K_n = ActiveVariable(self._reset_K_n,
"Coefficient K_n_0")
# Dependencies
self._xlambda.is_computed_by(
self.solid._rho_f, self.fluid._rho_f
......@@ -106,29 +105,21 @@ class BaseYosioka(BaseSphereFrequencyComposite):
)
self._A_n.is_computed_by(
self._xlambda, self.fluid._k_f,
self.solid._k_f, self._R_0, self._n
self.solid._k_f, self._R_0
)
self._B_n.is_computed_by(
self._xlambda, self.fluid._k_f,
self.solid._k_f, self._R_0, self._n
)
self._C_n.is_computed_by(
self._A_n, self._delta_n
)
self._D_n.is_computed_by(
self._B_n, self._delta_n
self.solid._k_f, self._R_0
)
self._delta_n.is_computed_by(
self._n, self._h, self.fluid._k_f
)
self._M_n_0.is_computed_by(
self._n, self._B_n, self._omega,
self._t, self.solid._k_f, self._R_0
)
self._K_n_0.is_computed_by(
self._n, self._B_n, self._omega,
self._t, self.solid._k_f, self._R_0
self.field._position, self.fluid._k_f
)
self._M_n.is_computed_by(self._omega,
self._t, self.solid._k_f, self._R_0
)
self._K_n.is_computed_by(self._omega,
self._t, self.solid._k_f, self._R_0
)
log.info(str(self))
log.debug(repr(self))
......@@ -186,34 +177,6 @@ class BaseYosioka(BaseSphereFrequencyComposite):
def F(self) -> float:
return self._F.value
@property
def A_n(self) -> complex:
return self._A_n.value
@property
def B_n(self) -> complex:
return self._B_n.value
@property
def C_n(self) -> complex:
return self._C_n.value
@property
def D_n(self) -> complex:
return self._D_n.value
@property
def delta_n(self) -> complex:
return self._delta_n.value
@property
def M_n_0(self) -> complex:
return self._M_n_0.value
@property
def K_n_0(self) -> complex:
return self._K_n_0.value
# -------------------------------------------------------------------------
# Wrappers for Independent Field Attributes
# -------------------------------------------------------------------------
......@@ -382,51 +345,118 @@ class BaseYosioka(BaseSphereFrequencyComposite):
# Dependent Variables Methods
# -------------------------------------------------------------------------
def _compute_A_n(self) -> complex:
@staticmethod
def _reset_A_n() -> list:
return []
@staticmethod
def _reset_B_n() -> list:
return []
@staticmethod
def _reset_C_n() -> list:
return []
@staticmethod
def _reset_D_n() -> list:
return []
@staticmethod
def _reset_delta_n() -> list:
return []
def A_n(self, n: int) -> complex:
if n < len(self._A_n.value):
return self._A_n.value[n]
else:
self._compute_A_n(n)
return self._A_n.value[n]
def B_n(self, n: int) -> complex:
if n < len(self._B_n.value):
return self._B_n.value[n]
else:
self._compute_B_n(n)
return self._B_n.value[n]
def C_n(self, n: int) -> complex:
if n < len(self._C_n.value):
return self._C_n.value[n]
else:
self._compute_C_n(n)
return self._C_n.value[n]
def D_n(self, n: int) -> complex:
if n < len(self._D_n.value):
return self._D_n.value[n]
else:
self._compute_D_n(n)
return self._D_n.value[n]
def delta_n(self, n: int) -> complex:
if n < len(self._delta_n.value):
return self._delta_n.value[n]
else:
self._compute_delta_n(n)
return self._delta_n.value[n]
def _compute_A_n(self, N: int):
"""
Compute A_n according to (22)
"""
s = self.k_s * self.R_0
f = self.k_f * self.R_0
term1 = self.xlambda * self.k_f * sp.jn(self.n, s) * sp.d_jn(self.n, f)
term2 = self.k_s * sp.d_jn(self.n, s) * sp.jn(self.n, f)
term3 = self.k_s * sp.d_jn(self.n, s) * sp.hn(self.n, f)
term4 = self.xlambda * self.k_f * sp.jn(self.n, s) * sp.d_hn(self.n, f)
out = (term1 - term2) / (term3 - term4)
return out
n_old = len(self._A_n.value)
for n in full_range(n_old, N):
s = self.k_s * self.R_0
f = self.k_f * self.R_0
term1 = self.xlambda * self.k_f * sp.jn(n, s) * sp.d_jn(n, f)
term2 = self.k_s * sp.d_jn(n, s) * sp.jn(n, f)
term3 = self.k_s * sp.d_jn(n, s) * sp.hn(n, f)
term4 = self.xlambda * self.k_f * sp.jn(n, s) * sp.d_hn(n, f)
out = (term1 - term2) / (term3 - term4)
self._A_n.value.append(out)
def _compute_B_n(self) -> complex:
def _compute_B_n(self, N: int):
"""
Compute B_n according to (23)
"""
s = self.k_s * self.R_0
f = self.k_f * self.R_0
out = 1j * self.k_f
out /= f ** 2
term1 = self.k_s * sp.d_jn(self.n, s) * sp.hn(self.n, f)
term2 = self.xlambda * self.k_f * sp.jn(self.n, s) * sp.d_hn(self.n, f)
out /= (term1 - term2)
return out
n_old = len(self._B_n.value)
for n in full_range(n_old, N):
s = self.k_s * self.R_0
f = self.k_f * self.R_0
out = 1j * self.k_f
out /= f ** 2
term1 = self.k_s * sp.d_jn(n, s) * sp.hn(n, f)
term2 = self.xlambda * self.k_f * sp.jn(n, s) * sp.d_hn(n, f)
out /= (term1 - term2)
self._B_n.value.append(out)
def _compute_C_n(self) -> complex:
def _compute_C_n(self, N: int):
"""
Compute C_n according to (31)
"""
return self._A_n.value * self.delta_n
n_old = len(self._C_n.value)
for n in full_range(n_old, N):
out = self.A_n(n) * self.delta_n(n)
self._C_n.value.append(out)
def _compute_D_n(self) -> complex:
def _compute_D_n(self, N: int):
"""
Compute D_n according to (32)
"""
return self._B_n.value * self.delta_n
n_old = len(self._D_n.value)
for n in full_range(n_old, N):
out = self.B_n(n) * self.delta_n(n)
self._D_n.value.append(out)
def _compute_delta_n(self) -> complex:
def _compute_delta_n(self, N: int):
"""
Compute delta_n according to (28)
"""
out = (-1) ** self.n * exp(1j * self.k_f * self.h)
out += exp(-1j * self.k_f * self.h)
return out
n_old = len(self._delta_n.value)
for n in full_range(n_old, N):
out = (-1) ** n * exp(1j * self.k_f * self.position)
out += exp(-1j * self.k_f * self.position)
self._delta_n.value.append(out)
def _compute_F(self) -> float:
"""
......@@ -468,21 +498,49 @@ class BaseYosioka(BaseSphereFrequencyComposite):
"""
return self.rho_s / self.rho_f
def _compute_M_n_0(self) -> complex:
@staticmethod
def _reset_K_n() -> list:
return []
@staticmethod
def _reset_M_n() -> list:
return []
def M_n(self, n: int) -> complex:
if n < len(self._M_n.value):
return self._M_n.value[n]
else:
self._compute_M_n(n)
return self._M_n.value[n]
def K_n(self, n: int) -> complex:
if n < len(self._K_n.value):
return self._K_n.value[n]
else:
self._compute_K_n(n)
return self._K_n.value[n]
def _compute_M_n(self, N: int):
"""
Compute M_n_0 according to (42)
Compute M_n according to (42)
"""
out = (-1j) ** self.n * self.B_n
out *= sp.jn(self.n, self.k_s * self.R_0)
return out
def _compute_K_n_0(self) -> complex:
n_old = len(self._M_n.value)
for n in full_range(n_old, N):
out = (-1j) ** n * self.B_n(n)
out *= sp.jn(n, self.k_s * self.R_0)
self._M_n.value.append(out)
def _compute_K_n(self, N: int):
"""
Compute K_n_0 according to (43)
Compute K_n according to (43)
"""
out = (-1j) ** self.n * self.B_n * self.k_s
out *= sp.d_jn(self.n, self.k_s * self.R_0)
return out
n_old = len(self._K_n.value)
for n in full_range(n_old, N):
out = (-1j) ** n * self.B_n(n) * self.k_s
out *= sp.d_jn(n, self.k_s * self.R_0)
self._K_n.value.append(out)
def _D(self, num: int, n: int) -> float:
"""Computes the coefficients D, where num takes values 1 or 2
......
......@@ -218,7 +218,7 @@ class ScatteringField(BaseYosioka, BaseScattering):
elif self.wave_type == "standing":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n
* self.delta_n
* self.delta_n(n)
* sp.jn(n, self.k_f * self.r))
else:
# No incident potential field
......@@ -244,7 +244,7 @@ class ScatteringField(BaseYosioka, BaseScattering):
elif self.wave_type == "standing":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n
* self.delta_n
* self.delta_n(n)
* self.k_f *
sp.d_jn(n, self.k_f * self.r))
else:
......@@ -271,7 +271,7 @@ class ScatteringField(BaseYosioka, BaseScattering):
elif self.wave_type == "standing":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n
* self.delta_n
* self.delta_n(n)
* sp.jn(n, self.k_f * self.r))
else:
# No incident potential field
......@@ -298,13 +298,13 @@ class ScatteringField(BaseYosioka, BaseScattering):
# (18)
if self.wave_type == "traveling":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n * self.A_n
* sp.hn(self.n, self.k_f * self.r))
(2 * n + 1) * (-1j) ** n * self.A_n(n)
* sp.hn(n, self.k_f * self.r))
# (29)
elif self.wave_type == "standing":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n * self.C_n
* sp.hn(self.n, self.k_f * self.r))
(2 * n + 1) * (-1j) ** n * self.C_n(n)
* sp.hn(n, self.k_f * self.r))
else:
# No incident potential field
coeff = 0
......@@ -322,15 +322,15 @@ class ScatteringField(BaseYosioka, BaseScattering):
# (18)
if self.wave_type == "traveling":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n * self.A_n
(2 * n + 1) * (-1j) ** n * self.A_n(n)
* self.k_f *
sp.d_hn(self.n, self.k_f * self.r))
sp.d_hn(n, self.k_f * self.r))
# (29)
elif self.wave_type == "standing":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n * self.C_n
(2 * n + 1) * (-1j) ** n * self.C_n(n)
* self.k_f *
sp.d_hn(self.n, self.k_f * self.r))
sp.d_hn(n, self.k_f * self.r))
else:
# No incident potential field
coeff = 0
......@@ -348,13 +348,13 @@ class ScatteringField(BaseYosioka, BaseScattering):
# (18)
if self.wave_type == "traveling":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n * self.A_n
* sp.hn(self.n, self.k_f * self.r))
(2 * n + 1) * (-1j) ** n * self.A_n(n)
* sp.hn(n, self.k_f * self.r))
# (29)
elif self.wave_type == "standing":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n * self.C_n
* sp.hn(self.n, self.k_f * self.r))
(2 * n + 1) * (-1j) ** n * self.C_n(n)
* sp.hn(n, self.k_f * self.r))
else:
# No incident potential field
coeff = 0
......@@ -379,13 +379,13 @@ class ScatteringField(BaseYosioka, BaseScattering):
# (19)
if self.wave_type == "traveling":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n * self.B_n
* sp.jn(self.n, self.k_s * self.r))
(2 * n + 1) * (-1j) ** n * self.B_n(n)
* sp.jn(n, self.k_s * self.r))
# (30)
elif self.wave_type == "standing":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n * self.D_n
* sp.jn(self.n, self.k_s * self.r))
(2 * n + 1) * (-1j) ** n * self.D_n(n)
* sp.jn(n, self.k_s * self.r))
else:
# No incident potential field
coeff = 0
......@@ -403,15 +403,15 @@ class ScatteringField(BaseYosioka, BaseScattering):
# (19)
if self.wave_type == "traveling":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n * self.B_n
(2 * n + 1) * (-1j) ** n * self.B_n(n)
* self.k_s *
sp.d_jn(self.n, self.k_s * self.r))
sp.d_jn(n, self.k_s * self.r))
# (30)
elif self.wave_type == "standing":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n * self.D_n
(2 * n + 1) * (-1j) ** n * self.D_n(n)
* self.k_s *
sp.d_jn(self.n, self.k_s * self.r))
sp.d_jn(n, self.k_s * self.r))
else:
# No incident potential field
coeff = 0
......@@ -430,13 +430,13 @@ class ScatteringField(BaseYosioka, BaseScattering):
# (19)
if self.wave_type == "traveling":