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

various bugfixes adressing issue #29 and #31

parent a8b5306b
Pipeline #88580 failed with stages
in 1 minute and 55 seconds
......@@ -2,7 +2,7 @@ from __future__ import annotations
from typing import Union
from gorkov import (
log, pi, full_range, exp,
log, pi, full_range,
BaseSphereFrequencyComposite,
Sphere, Frequency,
BackgroundField,
......@@ -62,10 +62,7 @@ class BaseYosioka(BaseSphereFrequencyComposite):
particle_type,
"Type of the particle")
# Independent variables:
self._theta = PassiveVariable(0, "angle of attack")
self._r = PassiveVariable(position, "Position")
self._omega = PassiveVariable(2 * pi * frequency, "Angular Frequency")
self._t = PassiveVariable(0, "time")
self._R_0 = PassiveVariable(radius, "Particle-Radius")
# Dependent variables
self._xlambda = ActiveVariable(self._compute_xlambda,
......@@ -105,12 +102,10 @@ class BaseYosioka(BaseSphereFrequencyComposite):
self._xlambda, self.fluid._k_f,
self.solid._k_f, self._R_0
)
self._M_n.is_computed_by(self._omega,
self._t, self.solid._k_f,
self._M_n.is_computed_by(self._omega, self.solid._k_f,
self._R_0, self._B_n
)
self._K_n.is_computed_by(self._omega,
self._t, self.solid._k_f,
self._K_n.is_computed_by(self._omega, self.solid._k_f,
self._R_0, self._B_n
)
......@@ -294,22 +289,6 @@ class BaseYosioka(BaseSphereFrequencyComposite):
# Wrappers for independent variables
# ----------------------------------
@property
def theta(self) -> float:
return self._theta.value
@theta.setter
def theta(self, value: float) -> None:
self._theta.value = value
@property
def r(self) -> float:
return self._r.value
@r.setter
def r(self, value: float) -> None:
self._r.value = value
@property
def R_0(self) -> float:
return self._R_0.value
......@@ -326,14 +305,6 @@ class BaseYosioka(BaseSphereFrequencyComposite):
def omega(self, value: float) -> None:
self._omega.value = value
@property
def t(self) -> float:
return self._t.value
@t.setter
def t(self, value: float) -> None:
self._t.value = value
# -------------------------------------------------------------------------
# Dependent Variables Methods
# -------------------------------------------------------------------------
......@@ -444,24 +415,20 @@ class BaseYosioka(BaseSphereFrequencyComposite):
"""
Compute M_n according to (42)
"""
n_old = len(self._M_n.value)
for n in full_range(n_old, N):
out = (-1j) ** n * self.B_n(n)
out = (1j) ** n * self.B_n(n)
out *= sp.jn(n, self.k_s * self.R_0)
out *= exp(1j * self.omega * self.t)
self._M_n.value.append(out)
def _compute_K_n(self, N: int) -> None:
"""
Compute K_n according to (43)
"""
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 = (1j) ** n * self.B_n(n) * self.k_s
out *= sp.d_jn(n, self.k_s * self.R_0)
out *= exp(1j * self.omega * self.t)
self._K_n.value.append(out)
def _D(self, num: int, n: int) -> float:
......@@ -495,21 +462,21 @@ class BaseYosioka(BaseSphereFrequencyComposite):
lambda function is given as an argument
"""
out = []
self.n = 0
n = 0
if self.n_max is None:
while True:
out.append(lambda_function(self.n))
if self.n >= 1:
error = abs((out[self.n] - out[self.n - 1]) /
out[self.n - 1])
out.append(lambda_function(n))
if n >= 1:
error = abs((out[n] - out[n - 1]) /
out[n - 1])
if error < threshold:
break
if self.n >= iterations:
if n >= iterations:
break
self.n += 1
n += 1
else:
for self.n in range(self.n_max):
out.append(lambda_function(self.n))
for n in range(self.n_max):
out.append(lambda_function(n))
return out
......
from typing import Optional, Union
from gorkov import (ActiveVariable,
Frequency, Sphere,
from gorkov import (Frequency, Sphere,
LegendreFunctions as leg,
SpecialFunctions as sp,
exp, log)
......@@ -24,160 +23,9 @@ class ScatteringField(BaseYosioka, BaseScattering):
p_0, position,
particle_type, wave_type, n_max)
# Dependent variables:
self._Phi_1 = ActiveVariable(self._compute_Phi_1,
"Velocity potential outside the sphere")
self._Phi_i = ActiveVariable(self._compute_Phi_i,
"Incident velocity potential")
self._Phi_s = ActiveVariable(self._compute_Phi_s,
"velocity potential of scattered waves")
self._Phi_star = ActiveVariable(self._compute_Phi_star,
"velocity potential inside particle")
self._d_Phi_1_dr = ActiveVariable(
self._compute_d_Phi_1_dr,
"Derivative of Phi_1 w.r.t. r")
self._d_Phi_1_dtheta = ActiveVariable(
self._compute_d_Phi_1_dtheta,
"Derivative of Phi_1 w.r.t. theta")
self._d_Phi_i_dr = ActiveVariable(
self._compute_d_Phi_i_dr,
"Derivative of Phi_i w.r.t. r")
self._d_Phi_i_dtheta = ActiveVariable(
self._compute_d_Phi_i_dtheta,
"Derivative of Phi_i w.r.t. theta")
self._d_Phi_s_dr = ActiveVariable(
self._compute_d_Phi_s_dr,
"Derivative of Phi_s w.r.t. r")
self._d_Phi_s_dtheta = ActiveVariable(
self._compute_d_Phi_s_dtheta,
"Derivative of Phi_s w.r.t. theta")
self._d_Phi_star_dr = ActiveVariable(
self._compute_d_Phi_star_dr,
"Derivative of Phi_star w.r.t. r")
self._d_Phi_star_dtheta = ActiveVariable(
self._compute_d_Phi_star_dtheta,
"Derivative of Phi_star w.r.t. theta")
self._radial_acoustic_fluid_velocity = ActiveVariable(
self.radial_acoustic_fluid_velocity, "radial acoustic velocity")
self._radial_particle_velocity = ActiveVariable(
self.radial_particle_velocity, "radial particle velocity")
self._tangential_acoustic_fluid_velocity = ActiveVariable(
self.tangential_acoustic_fluid_velocity,
"tangential acoustic velocity")
self._tangential_particle_velocity = ActiveVariable(
self.tangential_particle_velocity, "tangential particle velocity")
# Dependencies:
self._Phi_1.is_computed_by(
self._Phi_i, self._Phi_s)
self._d_Phi_1_dr.is_computed_by(
self._d_Phi_i_dr, self._d_Phi_s_dr)
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.fluid._k_f,
self._r, self._theta, self._t)
self._d_Phi_i_dtheta.is_computed_by(
self.fluid._k_f,
self._r, self._theta, self._t)
self._d_Phi_i_dr.is_computed_by(
self.fluid._k_f,
self._r, self._theta, self._t)
self._Phi_s.is_computed_by(
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.fluid._k_f,
self._r, self._theta, self._t)
self._d_Phi_s_dtheta.is_computed_by(
self._A_n, self.fluid._k_f,
self._r, self._theta, self._t)
self._Phi_star.is_computed_by(
self._B_n, self.fluid._k_f,
self._r, self._theta, self._t)
self._d_Phi_star_dr.is_computed_by(
self._B_n, self.fluid._k_f,
self._r, self._theta, self._t)
self._d_Phi_star_dtheta.is_computed_by(
self._B_n, self.fluid._k_f,
self._r, self._theta, self._t)
self._radial_particle_velocity.is_computed_by(
self._d_Phi_star_dr)
self._tangential_particle_velocity.is_computed_by(
self._d_Phi_star_dtheta, self._r)
self._radial_acoustic_fluid_velocity.is_computed_by(
self._d_Phi_1_dr)
self._tangential_acoustic_fluid_velocity.is_computed_by(
self._d_Phi_1_dtheta, self._r)
@property
def Phi_1(self) -> complex:
return self._Phi_1.value
@property
def d_Phi_1_dr(self) -> complex:
return self._d_Phi_s_dr.value
@property
def d_Phi_1_dtheta(self) -> complex:
return self._d_Phi_1_dtheta.value
@property
def Phi_i(self) -> complex:
return self._Phi_i.value
@property
def d_Phi_i_dr(self) -> complex:
return self._d_Phi_i_dr.value
@property
def d_Phi_i_dtheta(self) -> complex:
return self._d_Phi_i_dtheta.value
@property
def Phi_s(self) -> complex:
return self._Phi_s.value
@property
def d_Phi_s_dr(self) -> complex:
return self._d_Phi_s_dr.value
@property
def d_Phi_s_dtheta(self) -> complex:
return self._d_Phi_s_dtheta.value
@property
def Phi_star(self) -> complex:
return self._Phi_star.value
@property
def d_Phi_star_dr(self) -> complex:
return self._d_Phi_star_dr.value
@property
def d_Phi_star_dtheta(self) -> complex:
return self._d_Phi_star_dtheta.value
@property
def v_ac_r(self) -> float:
return self._radial_acoustic_fluid_velocity.value
@property
def v_p_r(self) -> float:
return self._radial_particle_velocity.value
@property
def v_ac_t(self) -> float:
return self._tangential_acoustic_fluid_velocity.value
@property
def v_p_t(self) -> float:
return self._tangential_particle_velocity.value
def __repr__(self):
return (
f'ScatteringField({self.f}, {self.r}, {self.rho_s}, '
f'ScatteringField({self.f}, {self.R_0}, {self.rho_s}, '
f'{self.rho_f}), {self.c_f}. {self.c_s}, {self.p_0}, '
f'{self.wave_type}, {self.particle_type}, '
f'{self.position}'
......@@ -187,21 +35,15 @@ class ScatteringField(BaseYosioka, BaseScattering):
# Velocity Potentials
# -------------------------------------------------------------------------
def _compute_Phi_1(self) -> complex:
def Phi_1(self, r: float, theta: float, t: float) -> complex:
"""
Compute velocity potential outside the sphere
according to (17)
"""
# (17)
return self.Phi_i + self.Phi_s
def _compute_d_Phi_1_dr(self) -> complex:
return self.d_Phi_i_dr + self.d_Phi_s_dr
return self.Phi_i(r, theta, t) + self.Phi_s(r, theta, t)
def _compute_d_Phi_1_dtheta(self) -> complex:
return self.d_Phi_i_dtheta + self.d_Phi_s_dtheta
def _compute_Phi_i(self) -> complex:
def Phi_i(self, r: float, theta: float, t: float) -> complex:
"""
Compute the incident velocity potential
according to (16) for travelling waves
......@@ -211,32 +53,13 @@ class ScatteringField(BaseYosioka, BaseScattering):
# (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.cos_poly(self.theta, coeff)
return out
def _compute_d_Phi_i_dr(self) -> complex:
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.cos_poly(self.theta, coeff)
return out
def _compute_d_Phi_i_dtheta(self) -> complex:
coeff = self.legendre_coeffs(lambda n:
self.field.A_in(n)
* sp.jn(n, self.k_f * self.r))
* sp.jn(n, self.k_f * r))
out = exp(-1j * (self.omega * self.t))
out *= leg.first_cos_poly(self.theta, coeff)
out = exp(-1j * (self.omega * t))
out *= leg.cos_poly(theta, coeff)
return out
def _compute_Phi_s(self) -> complex:
def Phi_s(self, r: float, theta: float, t: float) -> complex:
"""
Compute the scattered velocity potential
according to (18) for travelling waves
......@@ -246,32 +69,13 @@ class ScatteringField(BaseYosioka, BaseScattering):
# (18) & (29)
coeff = self.legendre_coeffs(lambda n:
self.A_n(n)
* sp.hn(n, self.k_f * self.r))
* sp.hn(n, self.k_f * r))
out = exp(-1j * (self.omega * self.t))
out *= leg.cos_poly(self.theta, coeff)
out = exp(-1j * (self.omega * t))
out *= leg.cos_poly(theta, coeff)
return out
def _compute_d_Phi_s_dr(self) -> complex:
coeff = self.legendre_coeffs(lambda 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.cos_poly(self.theta, coeff)
return out
def _compute_d_Phi_s_dtheta(self) -> complex:
coeff = self.legendre_coeffs(lambda 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_Phi_star(self) -> complex:
def Phi_star(self, r: float, theta: float, t: float) -> complex:
"""
Compute the velocity potential inside the sphere
according to (19) for travelling waves
......@@ -281,29 +85,10 @@ class ScatteringField(BaseYosioka, BaseScattering):
# (19) & (30)
coeff = self.legendre_coeffs(lambda n:
self.B_n(n)
* sp.jn(n, self.k_s * self.r))
* sp.jn(n, self.k_s * r))
out = exp(-1j * (self.omega * self.t))
out *= leg.cos_poly(self.theta, coeff)
return out
def _compute_d_Phi_star_dr(self) -> complex:
coeff = self.legendre_coeffs(lambda 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.cos_poly(self.theta, coeff)
return out
def _compute_d_Phi_star_dtheta(self) -> complex:
coeff = self.legendre_coeffs(lambda 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)
out = exp(-1j * (self.omega * t))
out *= leg.cos_poly(theta, coeff)
return out
# -------------------------------------------------------------------------
......@@ -321,10 +106,25 @@ class ScatteringField(BaseYosioka, BaseScattering):
:type t: float
:rtype: float
"""
self.t = t
self.theta = theta
self.r = r
return self.d_Phi_1_dr.real
coeff = self.legendre_coeffs(
lambda n:
self.field.A_in(n) *
self.k_f * sp.d_jn(n, self.k_f * r))
d_Phi_i_dr = leg.cos_poly(theta, coeff)
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 * r))
d_Phi_s_dr = leg.cos_poly(theta, coeff)
d_Phi_1_dr = (d_Phi_i_dr + d_Phi_s_dr)
d_Phi_1_dr *= exp(-1j * (self.omega * t))
return d_Phi_1_dr.real
def tangential_acoustic_fluid_velocity(self, r: float, theta: float,
t: float) -> float:
......@@ -337,10 +137,27 @@ class ScatteringField(BaseYosioka, BaseScattering):
:type t: float
:rtype: float
"""
self.t = t
self.theta = theta
self.r = r
return self.d_Phi_1_dtheta.real / self.r
coeff = self.legendre_coeffs(
lambda n:
self.field.A_in(n) *
sp.jn(n, self.k_f * r))
d_Phi_i_dtheta = leg.first_cos_poly(theta, coeff)
coeff = self.legendre_coeffs(
lambda n:
self.field.A_in(n) *
self.A_n(n) *
sp.hn(n, self.k_f * r))
d_Phi_s_dtheta = leg.first_cos_poly(theta, coeff)
d_Phi_1_dtheta = (d_Phi_i_dtheta + d_Phi_s_dtheta)
d_Phi_1_dtheta *= exp(-1j * (self.omega * t))
return d_Phi_1_dtheta.real / r
def radial_particle_velocity(self, r: float, theta: float,
t: float) -> float:
......@@ -353,10 +170,16 @@ class ScatteringField(BaseYosioka, BaseScattering):
:type t: float
:rtype: float
"""
self.t = t
self.theta = theta
self.r = r
return self.d_Phi_star_dr.real
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 * r))
d_Phi_star_dr = leg.cos_poly(theta, coeff)
d_Phi_star_dr *= exp(-1j * (self.omega * t))
return d_Phi_star_dr.real
def tangential_particle_velocity(self, r: float, theta: float,
t: float) -> float:
......@@ -372,21 +195,18 @@ class ScatteringField(BaseYosioka, BaseScattering):
:type t: float
:rtype: float
"""
self.t = t
self.theta = theta
self.r = r
return self.d_Phi_star_dtheta.real / self.r
def d_Phi1_dt(self, r: float, theta: float, t: float) -> float:
self.r = r
self.theta = theta
self.t = t
epsilon = t / 1000.0
out = self.Phi_1
self.t -= epsilon
out -= self.Phi_1
out /= epsilon
return out.real
coeff = self.legendre_coeffs(
lambda n:
self.field.A_in(n) *
self.B_n(n) *
sp.jn(n, self.k_s * r))
d_Phi_star_dtheta = leg.first_cos_poly(theta, coeff)
d_Phi_star_dtheta *= exp(-1j * (self.omega * t))
return d_Phi_star_dtheta.real / r
if __name__ == '__main__':
......
import unittest
import numpy as np
from gorkov.solutions import King1934, Gorkov1962
from math import isclose
from tests.basetest_composite import BaseTest
from gorkov import pi, sin, integrate
from gorkov import Yosioka1955
from gorkov import (
pi, sin, integrate,
Yosioka1955, King1934, Gorkov1962
)
class TestYosiokaARF(BaseTest):
def assertIsClose(self, a: float, b: float) -> None:
rel_tolerance = 1e-6
self.assertTrue(isclose(a, b, rel_tol=rel_tolerance),
msg=f'expected: {b}, got: {a}')
def assertIsNotClose(self, a: float, b: float) -> None:
rel_tolerance = 1e-6
self.assertFalse(isclose(a, b, rel_tol=rel_tolerance),
msg=f'expected: {b}, got: {a}')
def setUp(self) -> None:
# Geometry
......@@ -145,6 +156,28 @@ class TestYosiokaARF(BaseTest):
self.cls.t = 5e-6
print(integrate(self.cls.time_averaging_uv, 0, 1.0 / self.f))
def test_compute_coefficients(self) -> None:
"""
check if coefficients are computed correctly
"""
real_coeff = dict({self.cls.A_n(0).real: -2.1259984038165364e-14,
self.cls.B_n(0).real: -0.9523743230446932,
self.cls.K_n(0).real: 2.7849320017549104,
self.cls.M_n(0).real: -0.9523046982176053