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

Improving analytical derivative

parent 98b64a49
Pipeline #88019 canceled with stages
in 1 minute and 55 seconds
...@@ -278,11 +278,15 @@ class ScatteringField(BaseYosioka, BaseScattering): ...@@ -278,11 +278,15 @@ class ScatteringField(BaseYosioka, BaseScattering):
coeff = 0 coeff = 0
raise NotImplementedError raise NotImplementedError
out = exp(1j * (self.omega * self.t)) out = exp(1j * (self.omega * self.t))
degrees = len(coeff)
out *= sum([coeff[n] * (lpmv(2, n, cos(self.theta))) p = [lpmv(1, n, cos(self.theta))
for n in range(degrees)]) for n in range(len(coeff))]
out *= -sin(self.theta) p_der = []
out /= self.r for i, val in enumerate(p):
if i > 0:
p_der.append(i * val * coeff[i - 1])
out *= sum(p_der)
out *= -sin(self.theta) / self.r
return out return out
def _compute_Phi_s(self) -> complex: def _compute_Phi_s(self) -> complex:
...@@ -356,11 +360,14 @@ class ScatteringField(BaseYosioka, BaseScattering): ...@@ -356,11 +360,14 @@ class ScatteringField(BaseYosioka, BaseScattering):
coeff = 0 coeff = 0
raise NotImplementedError raise NotImplementedError
out = exp(1j * (self.omega * self.t)) out = exp(1j * (self.omega * self.t))
degrees = len(coeff) p = [lpmv(1, n, cos(self.theta))
out *= sum([coeff[n] * (lpmv(2, n, cos(self.theta))) for n in range(len(coeff))]
for n in range(degrees)]) p_der = []
out *= -sin(self.theta) for i, val in enumerate(p):
out /= self.r if i > 0:
p_der.append(i * val * coeff[i - 1])
out *= sum(p_der)
out *= -sin(self.theta) / self.r
return out return out
def _compute_Phi_star(self) -> complex: def _compute_Phi_star(self) -> complex:
...@@ -419,6 +426,7 @@ class ScatteringField(BaseYosioka, BaseScattering): ...@@ -419,6 +426,7 @@ class ScatteringField(BaseYosioka, BaseScattering):
according to (19) for travelling waves according to (19) for travelling waves
and (30) for standing waves and (30) for standing waves
""" """
# (19) # (19)
if self.wave_type == "traveling": if self.wave_type == "traveling":
coeff = self.legendre_coeffs(lambda n: coeff = self.legendre_coeffs(lambda n:
...@@ -433,13 +441,15 @@ class ScatteringField(BaseYosioka, BaseScattering): ...@@ -433,13 +441,15 @@ class ScatteringField(BaseYosioka, BaseScattering):
# No incident potential field # No incident potential field
coeff = 0 coeff = 0
raise NotImplementedError raise NotImplementedError
out = exp(1j * (self.omega * self.t)) out = exp(1j * (self.omega * self.t))
degrees = len(coeff) p = [lpmv(1, n, cos(self.theta))
out *= sum([coeff[n] * (lpmv(2, n, cos(self.theta))) for n in range(len(coeff))]
for n in range(degrees)]) p_der = []
out *= -sin(self.theta) for i, val in enumerate(p):
out /= self.r if i > 0:
p_der.append(i * val * coeff[i - 1])
out *= sum(p_der)
out *= -sin(self.theta) / self.r
return out return out
# ------------------------------------------------------------------------- # -------------------------------------------------------------------------
......
...@@ -37,7 +37,7 @@ class TestYosiokaARF(BaseTest): ...@@ -37,7 +37,7 @@ class TestYosiokaARF(BaseTest):
self.wave_type = 'traveling' self.wave_type = 'traveling'
self.inf_type = 'radius' self.inf_type = 'radius'
self.pos = 1e-4 self.pos = 1e-4
n_max = 3 n_max = 1e10
self.cls = Yosioka1955.ARF(frequency=self.f, self.cls = Yosioka1955.ARF(frequency=self.f,
radius=self.R_0, rho_s=self.rho_s, radius=self.R_0, rho_s=self.rho_s,
...@@ -85,13 +85,14 @@ class TestYosiokaARF(BaseTest): ...@@ -85,13 +85,14 @@ class TestYosiokaARF(BaseTest):
'traveling', self.p0, self.pos) 'traveling', self.p0, self.pos)
self.cls_list = [ self.cls_list = [
self.cls, self.cls
self.cls_num, # ,
self.cls_special, # self.cls_num,
self.cls_bubble, # self.cls_special,
self.cls_sphere, # self.cls_bubble,
self.clsKing, # self.cls_sphere,
self.clsGorkov # self.clsKing,
# self.clsGorkov
] ]
def randomly_change_parameters(self) -> None: def randomly_change_parameters(self) -> None:
...@@ -148,8 +149,8 @@ class TestYosiokaARF(BaseTest): ...@@ -148,8 +149,8 @@ class TestYosiokaARF(BaseTest):
def test_arf(self) -> None: def test_arf(self) -> None:
n = 0 n = 0
while (n < 5): while n < 1:
self.change_and_assign() # self.change_and_assign()
for cls in self.cls_list: for cls in self.cls_list:
print(cls.acoustic_radiation_force()) print(cls.acoustic_radiation_force())
print("") print("")
......
...@@ -96,7 +96,7 @@ class TestScatteringField(unittest.TestCase): ...@@ -96,7 +96,7 @@ class TestScatteringField(unittest.TestCase):
self.cls.wave_type = 'traveling' self.cls.wave_type = 'traveling'
self.cls.n_max = 3 self.cls.n_max = 3
self.cls.r = 1e-6 self.cls.r = 1e-6
acc = 20 acc = 3
theta = np.linspace(0, 2 * pi, acc) theta = np.linspace(0, 2 * pi, acc)
t = np.linspace(0, 5e-6, acc) t = np.linspace(0, 5e-6, acc)
phi = np.zeros((len(theta), len(t))) phi = np.zeros((len(theta), len(t)))
...@@ -112,21 +112,22 @@ class TestScatteringField(unittest.TestCase): ...@@ -112,21 +112,22 @@ class TestScatteringField(unittest.TestCase):
def test_velocities(self) -> None: def test_velocities(self) -> None:
self.cls.wave_type = 'traveling' self.cls.wave_type = 'traveling'
self.cls.n_max = 10 self.cls.n_max = 5
self.cls.t = 1e-6 self.cls.t = 0
acc = 10 acc = 6
theta = np.linspace(0, 2 * pi, acc) theta = np.linspace(0, 2 * pi, acc)
r = np.linspace(1e-9, 1e-6, acc) r = np.linspace(1e-9, 1e-6, acc)
vel = np.zeros((len(theta), len(r))) vel = np.zeros([acc, acc])
for ix, dx in enumerate(theta): for i_t, val_t in enumerate(theta):
self.cls.theta = dx self.theta = val_t
for ir, dr in enumerate(r): for i_r, val_r in enumerate(r):
self.cls.r = dr self.r = val_r
vel[ix][ir] = self.cls.tangential_particle_velocity( vel[i_r][i_t] = self.cls.tangential_particle_velocity(
self.cls.r, self.cls.theta, self.cls.t) r=self.r, theta=self.theta, t=self.cls.t)
theta, r = np.meshgrid(theta, r) print(vel)
fig, ax = plt.subplots(subplot_kw={"projection": "3d"}) fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
ax.plot_surface(theta, r, vel, antialiased=True) theta, r = np.meshgrid(theta, r)
ax.plot_surface(theta, r, vel)
plt.show() plt.show()
......
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