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 76ccac01 authored by Cyrill Mast's avatar Cyrill Mast
Browse files

Bugfix: issue #29 (Potentials)

parent 44c6df98
......@@ -107,7 +107,7 @@ class ARF(BaseYosioka, BaseARF, BaseSphereFrequencyComposite):
while n < self.n_max:
# (45)
s1 = (self.K_n(n) * conj(self.K_n(n + 1))).real
s1 *= self.n + 1
s1 *= n + 1
out[0] += s1
# (46)
......
......@@ -2,7 +2,7 @@ from __future__ import annotations
from typing import Union
from gorkov import (
log, pi, exp, full_range,
log, pi, full_range,
BaseSphereFrequencyComposite,
Sphere, Frequency,
BackgroundField,
......@@ -57,14 +57,12 @@ class BaseYosioka(BaseSphereFrequencyComposite):
# Independent variables
self.n_max = n_max
self.n = 0
self._particle_type = PassiveVariable(
particle_type,
"Type of the particle")
# Independent variables:
self._theta = PassiveVariable(0, "angle of attack")
self._n = PassiveVariable(self.n, "Index of iteration")
self._r = PassiveVariable(position, "Position")
self._omega = PassiveVariable(2 * pi * frequency, "Angular Frequency")
self._t = PassiveVariable(0, "time")
......@@ -77,19 +75,15 @@ class BaseYosioka(BaseSphereFrequencyComposite):
self._F = ActiveVariable(self._compute_F,
"density-compressibility factor")
self._A_n = ActiveVariable(self._reset_A_n,
"Coefficient A_n for traveling field")
"Coefficient A_n"
" for scattered potential")
self._B_n = ActiveVariable(self._reset_B_n,
"Coefficient B_n for traveling field")
self._C_n = ActiveVariable(self._reset_C_n,
"Coefficient C_n for standing field")
self._D_n = ActiveVariable(self._reset_D_n,
"Coefficient C_n for standing field")
self._delta_n = ActiveVariable(self._reset_delta_n,
"Coefficient delta_n for calculations")
"Coefficient B_n"
" for potential inside particle")
self._M_n = ActiveVariable(self._reset_M_n,
"Coefficient M_n_0")
"Coefficient M_n")
self._K_n = ActiveVariable(self._reset_K_n,
"Coefficient K_n_0")
"Coefficient K_n")
# Dependencies
self._xlambda.is_computed_by(
self.solid._rho_f, self.fluid._rho_f
......@@ -111,9 +105,6 @@ class BaseYosioka(BaseSphereFrequencyComposite):
self._xlambda, self.fluid._k_f,
self.solid._k_f, self._R_0
)
self._delta_n.is_computed_by(
self.field._position, self.fluid._k_f
)
self._M_n.is_computed_by(self._omega,
self._t, self.solid._k_f, self._R_0
)
......@@ -353,18 +344,6 @@ class BaseYosioka(BaseSphereFrequencyComposite):
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]
......@@ -379,27 +358,6 @@ class BaseYosioka(BaseSphereFrequencyComposite):
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)
......@@ -430,34 +388,6 @@ class BaseYosioka(BaseSphereFrequencyComposite):
out /= (term1 - term2)
self._B_n.value.append(out)
def _compute_C_n(self, N: int):
"""
Compute C_n according to (31)
"""
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, N: int):
"""
Compute D_n according to (32)
"""
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, N: int):
"""
Compute delta_n according to (28)
"""
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:
"""
Compute F according to (60), (63), (75)
......
......@@ -77,31 +77,31 @@ class ScatteringField(BaseYosioka, BaseScattering):
self._d_Phi_1_dtheta.is_computed_by(
self._d_Phi_i_dtheta, self._d_Phi_s_dtheta)
self._Phi_i.is_computed_by(
self._delta_n, self.fluid._k_f,
self.fluid._k_f,
self._r, self._theta, self._t)
self._d_Phi_i_dtheta.is_computed_by(
self._delta_n, self.fluid._k_f,
self.fluid._k_f,
self._r, self._theta, self._t)
self._d_Phi_i_dr.is_computed_by(
self._delta_n, self.fluid._k_f,
self.fluid._k_f,
self._r, self._theta, self._t)
self._Phi_s.is_computed_by(
self._A_n, self._C_n, self.fluid._k_f,
self._A_n, self.fluid._k_f,
self._r, self._theta, self._t)
self._d_Phi_s_dr.is_computed_by(
self._A_n, self._C_n, self.fluid._k_f,
self._A_n, self.fluid._k_f,
self._r, self._theta, self._t)
self._d_Phi_s_dtheta.is_computed_by(
self._A_n, self._C_n, self.fluid._k_f,
self._A_n, self.fluid._k_f,
self._r, self._theta, self._t)
self._Phi_star.is_computed_by(
self._A_n, self._D_n, self.fluid._k_f,
self._B_n, self.fluid._k_f,
self._r, self._theta, self._t)
self._d_Phi_star_dr.is_computed_by(
self._A_n, self._D_n, self.fluid._k_f,
self._B_n, self.fluid._k_f,
self._r, self._theta, self._t)
self._d_Phi_star_dtheta.is_computed_by(
self._A_n, self._D_n, self.fluid._k_f,
self._B_n, self.fluid._k_f,
self._r, self._theta, self._t)
self._radial_particle_velocity.is_computed_by(
......@@ -207,76 +207,35 @@ class ScatteringField(BaseYosioka, BaseScattering):
"""
Compute the incident velocity potential
according to (16) for travelling waves
and (27) for standing waves
and (27) for standing waves.
Amplitude of the potential is multiplied.
"""
# (16)
if self.wave_type == "traveling":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n
* sp.jn(n, self.k_f * self.r))
# (27)
elif self.wave_type == "standing":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n
* self.delta_n(n)
* sp.jn(n, self.k_f * self.r))
else:
# No incident potential field
coeff = 0
# (16) & (27)
coeff = self.legendre_coeffs(lambda n:
self.field.A_in(n)
* sp.jn(n, self.k_f * self.r))
out = exp(1j * (self.omega * self.t))
out *= leg.first_cos_poly(self.theta,
coeff)
return out
def _compute_d_Phi_i_dr(self) -> complex:
"""
Compute the incident velocity potential
according to (16) for travelling waves
and (27) for standing waves
"""
# (16)
if self.wave_type == "traveling":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n
* self.k_f *
sp.d_jn(n, self.k_f * self.r))
# (27)
elif self.wave_type == "standing":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n
* self.delta_n(n)
* self.k_f *
sp.d_jn(n, self.k_f * self.r))
else:
# No incident potential field
coeff = 0
raise NotImplementedError
coeff = self.legendre_coeffs(lambda n:
self.field.A_in(n)
* self.k_f *
sp.d_jn(n, self.k_f * self.r))
out = exp(1j * (self.omega * self.t))
out *= leg.first_cos_poly(self.theta,
coeff)
return out
def _compute_d_Phi_i_dtheta(self) -> complex:
"""
Compute the incident velocity potential
according to (16) for travelling waves
and (27) for standing waves
"""
# (16)
if self.wave_type == "traveling":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n
* sp.jn(n, self.k_f * self.r))
# (27)
elif self.wave_type == "standing":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n
* self.delta_n(n)
* sp.jn(n, self.k_f * self.r))
else:
# No incident potential field
coeff = 0
raise NotImplementedError
coeff = self.legendre_coeffs(lambda n:
self.field.A_in(n)
* sp.jn(n, self.k_f * self.r))
out = exp(1j * (self.omega * self.t))
p = [lpmv(1, n, cos(self.theta))
......@@ -294,71 +253,35 @@ class ScatteringField(BaseYosioka, BaseScattering):
Compute the scattered velocity potential
according to (18) for travelling waves
and (29) for standing waves
Amplitude of the potential is multiplied.
"""
# (18)
if self.wave_type == "traveling":
coeff = self.legendre_coeffs(lambda n:
(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(n)
* sp.hn(n, self.k_f * self.r))
else:
# No incident potential field
coeff = 0
raise NotImplementedError
# (18) & (29)
coeff = self.legendre_coeffs(lambda n:
self.field.A_in(n)
* self.A_n(n)
* sp.hn(n, self.k_f * self.r))
out = exp(1j * (self.omega * self.t))
out *= leg.first_cos_poly(self.theta, coeff)
return out
def _compute_d_Phi_s_dr(self) -> complex:
"""
Compute the scattered velocity potential
according to (18) for travelling waves
and (29) for standing waves
"""
# (18)
if self.wave_type == "traveling":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n * self.A_n(n)
* self.k_f *
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(n)
* self.k_f *
sp.d_hn(n, self.k_f * self.r))
else:
# No incident potential field
coeff = 0
raise NotImplementedError
coeff = self.legendre_coeffs(lambda n:
self.field.A_in(n) *
self.A_n(n)
* self.k_f *
sp.d_hn(n, self.k_f * self.r))
out = exp(1j * (self.omega * self.t))
out *= leg.first_cos_poly(self.theta, coeff)
return out
def _compute_d_Phi_s_dtheta(self) -> complex:
"""
Compute the scattered velocity potential
according to (18) for travelling waves
and (29) for standing waves
"""
# (18)
if self.wave_type == "traveling":
coeff = self.legendre_coeffs(lambda n:
(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(n)
* sp.hn(n, self.k_f * self.r))
else:
# No incident potential field
coeff = 0
raise NotImplementedError
coeff = self.legendre_coeffs(lambda n:
self.field.A_in(n) *
self.A_n(n)
* sp.hn(n, self.k_f * self.r))
out = exp(1j * (self.omega * self.t))
p = [lpmv(1, n, cos(self.theta))
for n in range(len(coeff))]
......@@ -375,72 +298,37 @@ class ScatteringField(BaseYosioka, BaseScattering):
Compute the velocity potential inside the sphere
according to (19) for travelling waves
and (30) for standing waves
Amplitude of the potential is multiplied.
"""
# (19)
if self.wave_type == "traveling":
coeff = self.legendre_coeffs(lambda n:
(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(n)
* sp.jn(n, self.k_s * self.r))
else:
# No incident potential field
coeff = 0
raise NotImplementedError
# (19) & (30)
coeff = self.legendre_coeffs(lambda n:
self.field.A_in(n) *
self.B_n(n)
* sp.jn(n, self.k_s * self.r))
out = exp(1j * (self.omega * self.t))
out *= leg.first_cos_poly(self.theta, coeff)
return out
def _compute_d_Phi_star_dr(self) -> complex:
"""
Compute the velocity potential inside the sphere
according to (19) for travelling waves
and (30) for standing waves
"""
# (19)
if self.wave_type == "traveling":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n * self.B_n(n)
* self.k_s *
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(n)
* self.k_s *
sp.d_jn(n, self.k_s * self.r))
else:
# No incident potential field
coeff = 0
raise NotImplementedError
coeff = self.legendre_coeffs(lambda n:
self.field.A_in(n)
* self.B_n(n)
* self.k_s *
sp.d_jn(n, self.k_s * self.r))
out = exp(1j * (self.omega * self.t))
out *= leg.first_cos_poly(self.theta, coeff)
return out
def _compute_d_Phi_star_dtheta(self) -> complex:
"""
Compute the velocity potential inside the sphere
according to (19) for travelling waves
and (30) for standing waves
"""
# (19)
if self.wave_type == "traveling":
coeff = self.legendre_coeffs(lambda n:
(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(n)
* sp.jn(n, self.k_s * self.r))
else:
# No incident potential field
coeff = 0
raise NotImplementedError
coeff = self.legendre_coeffs(lambda n:
self.field.A_in(n)
* self.B_n(n)
* sp.jn(n, self.k_s * self.r))
out = exp(1j * (self.omega * self.t))
p = [lpmv(1, n, cos(self.theta))
for n in range(len(coeff))]
......
......@@ -52,7 +52,7 @@ class TestScatteringField(unittest.TestCase):
self.pos = 1e-5
self.theta = pi
self.n_max = 3
self.n_max = 5
self.cls = ScatteringField(self.f, self.R_0,
self.rho_s, self.rho_f, self.c_f,
......@@ -96,7 +96,7 @@ class TestScatteringField(unittest.TestCase):
self.cls.wave_type = 'traveling'
self.cls.n_max = 3
self.cls.r = 1e-6
acc = 3
acc = 20
theta = np.linspace(0, 2 * pi, acc)
t = np.linspace(0, 5e-6, acc)
phi = np.zeros((len(theta), len(t)))
......@@ -114,7 +114,7 @@ class TestScatteringField(unittest.TestCase):
self.cls.wave_type = 'traveling'
self.cls.n_max = 5
self.cls.t = 0
acc = 6
acc = 20
theta = np.linspace(0, 2 * pi, acc)
r = np.linspace(1e-9, 1e-6, acc)
vel = np.zeros([acc, acc])
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment