Skip to content

Commit 374694a

Browse files
MasashiSodeMasashiSode
committed
Merge remote-tracking branch 'origin/develop' into develop
2 parents 670a833 + 7d7fd6c commit 374694a

7 files changed

Lines changed: 854 additions & 15 deletions

File tree

‎examples/simple_plane_wave/simple_plane_wave_velocity_source.py‎

Lines changed: 197 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,197 @@
1+
"""Simple plane wave transmit example using a velocity (u0) source.
2+
3+
A velocity source drives the particle-velocity equation directly, as opposed to
4+
the hard pressure source (p0) which drives the pressure equation. For a plane
5+
wave propagating in the depth (x) direction the relevant component is u0.
6+
7+
The acoustic impedance relation p = rho · c · u links the two formulations, so
8+
the velocity amplitude is scaled as
9+
10+
u_amp = p_amp / (ρ₀ · c₀)
11+
12+
Everything else (grid, medium, sensor, solver) is identical to
13+
simple_plane_wave.py so the two examples can be run side-by-side for comparison.
14+
"""
15+
16+
import logging
17+
from pathlib import Path
18+
19+
import numpy as np
20+
21+
import fullwave
22+
from fullwave.utils import plot_utils, signal_process
23+
from fullwave.utils.coordinates import map_to_coords
24+
25+
26+
def main() -> None:
27+
"""Run simple plane wave transmit example with a velocity source."""
28+
logging.getLogger("__main__").setLevel(logging.INFO)
29+
30+
#
31+
# --- working directory ---
32+
#
33+
work_dir = Path("./outputs/") / "simple_plane_wave_velocity_source"
34+
work_dir.mkdir(parents=True, exist_ok=True)
35+
36+
#
37+
# --- computational grid ---
38+
#
39+
domain_size = (3e-2, 2e-2) # (depth, lateral) in metres
40+
f0 = 3e6 # centre frequency [Hz]
41+
c0 = 1540.0 # background sound speed [m/s]
42+
rho0 = 1000.0 # background density [kg/m³]
43+
duration = domain_size[0] / c0 * 2 # two-way travel time across depth
44+
grid = fullwave.Grid(
45+
domain_size=domain_size,
46+
f0=f0,
47+
duration=duration,
48+
c0=c0,
49+
)
50+
51+
#
52+
# --- acoustic medium ---
53+
#
54+
sound_speed_map = c0 * np.ones((grid.nx, grid.ny)) # m/s
55+
density_map = rho0 * np.ones((grid.nx, grid.ny)) # kg/m³
56+
alpha_coeff_map = 0.5 * np.ones((grid.nx, grid.ny)) # dB/(MHz^y cm)
57+
alpha_power_map = 1.0 * np.ones((grid.nx, grid.ny)) # power-law exponent
58+
beta_map = 0.0 * np.ones((grid.nx, grid.ny)) # nonlinearity coefficient
59+
60+
# embed a scatterer with different acoustic properties
61+
obj_x_start = grid.nx // 3
62+
obj_x_end = 2 * grid.nx // 3
63+
obj_y_start = grid.ny // 3
64+
obj_y_end = 2 * grid.ny // 3
65+
66+
sound_speed_map[obj_x_start:obj_x_end, obj_y_start:obj_y_end] = 1600
67+
density_map[obj_x_start:obj_x_end, obj_y_start:obj_y_end] = 1100
68+
alpha_coeff_map[obj_x_start:obj_x_end, obj_y_start:obj_y_end] = 0.75
69+
alpha_power_map[obj_x_start:obj_x_end, obj_y_start:obj_y_end] = 1.1
70+
71+
medium = fullwave.Medium(
72+
grid=grid,
73+
sound_speed=sound_speed_map,
74+
density=density_map,
75+
alpha_coeff=alpha_coeff_map,
76+
alpha_power=alpha_power_map,
77+
beta=beta_map,
78+
)
79+
medium.plot(export_path=work_dir / "medium.png")
80+
81+
#
82+
# --- velocity source ---
83+
#
84+
# Use the u-component (depth / x direction) to drive a downward-travelling
85+
# plane wave. The source occupies the top `element_thickness_px` rows of the
86+
# grid, matching the pressure-source layout in simple_plane_wave.py.
87+
#
88+
# Velocity amplitude is derived from the target pressure amplitude via the
89+
# plane-wave impedance relation: u_amp = p_amp / (ρ₀ · c₀)
90+
#
91+
p_amp = 1e5 # target pressure amplitude [Pa]
92+
u_amp = p_amp / (rho0 * c0) # corresponding velocity amplitude [m/s]
93+
94+
ncycles = 2
95+
drop_off = 2
96+
# element_thickness_px = 3
97+
98+
# Build the coordinate array for the velocity source layer
99+
# small velocity source at the center of the domain
100+
101+
source_width_px_x = 2
102+
source_width_px_y = 2
103+
u_mask = np.zeros((grid.nx, grid.ny), dtype=bool)
104+
u_mask[
105+
grid.nx // 2 - source_width_px_x // 2 : grid.nx // 2 + source_width_px_x // 2,
106+
grid.ny // 2 - source_width_px_y // 2 : grid.ny // 2 + source_width_px_y // 2,
107+
] = True
108+
109+
coords_u = map_to_coords(u_mask) # shape [n_sources_u, 2]
110+
111+
# Build the u0 signal matrix [n_sources_u, nt]
112+
u0 = np.zeros((coords_u.shape[0], grid.nt))
113+
114+
u0_vec = fullwave.utils.pulse.gaussian_modulated_sinusoidal_signal(
115+
nt=grid.nt,
116+
f0=f0,
117+
duration=duration,
118+
ncycles=ncycles,
119+
drop_off=drop_off,
120+
p0=u_amp, # amplitude in m/s
121+
)
122+
u0[:, :] = u0_vec
123+
124+
# for i_layer in range(element_thickness_px):
125+
# u0_vec = fullwave.utils.pulse.gaussian_modulated_sinusoidal_signal(
126+
# nt=grid.nt,
127+
# f0=f0,
128+
# duration=duration,
129+
# ncycles=ncycles,
130+
# drop_off=drop_off,
131+
# p0=u_amp, # amplitude in m/s
132+
# i_layer=i_layer,
133+
# dt_for_layer_delay=grid.dt,
134+
# cfl_for_layer_delay=grid.cfl,
135+
# )
136+
# n_y = coords_u.shape[0] // element_thickness_px
137+
# u0[n_y * i_layer : n_y * (i_layer + 1), :] = u0_vec
138+
139+
source = fullwave.Source(
140+
grid_shape=grid.shape,
141+
u0=u0,
142+
coords_u=coords_u,
143+
)
144+
145+
#
146+
# --- sensor ---
147+
#
148+
sensor_mask = np.ones((grid.nx, grid.ny), dtype=bool)
149+
sensor = fullwave.Sensor(mask=sensor_mask, sampling_modulus_time=7)
150+
151+
#
152+
# --- solver ---
153+
#
154+
fw_solver = fullwave.Solver(
155+
work_dir=work_dir,
156+
grid=grid,
157+
medium=medium,
158+
source=source,
159+
sensor=sensor,
160+
run_on_memory=False,
161+
use_exponential_attenuation=True,
162+
save_gpu_memory=True,
163+
path_fullwave_simulation_bin=Path(
164+
"/home/msode/workspace/lab_repos/fullwave-python-public/debug_solver_bin/fullwave2_2d_exponential_attenuation_multi_gpu",
165+
),
166+
)
167+
sensor_output = fw_solver.run()
168+
169+
#
170+
# --- visualisation ---
171+
#
172+
propagation_map = signal_process.reshape_whole_sensor_to_nt_nx_ny(
173+
sensor_output,
174+
grid,
175+
)
176+
p_max_plot = np.abs(propagation_map).max().item() / 4
177+
time_step = propagation_map.shape[0] // 3
178+
plot_utils.plot_array(
179+
propagation_map[time_step, :, :],
180+
aspect=propagation_map.shape[2] / propagation_map.shape[1],
181+
export_path=work_dir / "wave_propagation_snapshot.png",
182+
vmax=p_max_plot,
183+
vmin=-p_max_plot,
184+
)
185+
plot_utils.plot_wave_propagation_with_map(
186+
propagation_map=propagation_map,
187+
c_map=medium.sound_speed,
188+
rho_map=medium.density,
189+
export_name=work_dir / "wave_propagation_animation.mp4",
190+
vmax=p_max_plot,
191+
vmin=-p_max_plot,
192+
figsize=(4, 6),
193+
)
194+
195+
196+
if __name__ == "__main__":
197+
main()

‎fullwave/solver/binary_manager.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,7 @@
2828

2929
# Pinned release tag for the solver binaries.
3030
# Update this only when new binaries are uploaded to a GitHub release.
31-
BINARY_RELEASE_TAG = "fullwave_bin_v1.1"
31+
BINARY_RELEASE_TAG = "fullwave_bin_v1.2"
3232

3333

3434
def _download_url(filename: str, tag: str) -> str:

‎fullwave/solver/input_file_writer.py‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -169,6 +169,19 @@ def run(
169169
)
170170
if incoords_add is not None:
171171
self._queue_coords_write(simulation_dir / "icc_add.dat", incoords_add)
172+
for _vel_suffix, _vel_attr in (("u", "u0"), ("v", "v0"), ("w", "w0")):
173+
_signal_vel = getattr(self.source, _vel_attr, None)
174+
_incoords_vel = getattr(self.source, f"incoords_{_vel_suffix}", None)
175+
if _signal_vel is not None:
176+
self._queue_ic_write(
177+
simulation_dir / f"icmat_{_vel_suffix}.dat",
178+
np.transpose(_signal_vel),
179+
)
180+
if _incoords_vel is not None:
181+
self._queue_coords_write(
182+
simulation_dir / f"icc_{_vel_suffix}.dat",
183+
_incoords_vel,
184+
)
172185
self._queue_coords_write(simulation_dir / "icc.dat", self.source.incoords)
173186
self._copy_simulation_bin_file(simulation_dir)
174187

@@ -1047,6 +1060,10 @@ def _save_coords_params(self, simulation_dir: Path) -> None:
10471060
n_sources_add = getattr(self.source, "n_sources_add", 0)
10481061
if n_sources_add > 0:
10491062
var_list.append(("ncoords_add", n_sources_add))
1063+
for _vel_suffix in ("u", "v", "w"):
1064+
_n_vel = getattr(self.source, f"n_sources_{_vel_suffix}", 0)
1065+
if _n_vel > 0:
1066+
var_list.append((f"ncoords_{_vel_suffix}", _n_vel))
10501067
if self.is_3d:
10511068
var_list.extend(
10521069
[

‎fullwave/solver/pml_builder.py‎

Lines changed: 43 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -322,12 +322,33 @@ def __init__( # noqa: PLR0915
322322
if getattr(self.source_org, "incoords_add", None) is not None
323323
else None
324324
)
325+
incoords_u_ext = (
326+
self.source_org.incoords_u + self.num_boundary_points
327+
if getattr(self.source_org, "incoords_u", None) is not None
328+
else None
329+
)
330+
incoords_v_ext = (
331+
self.source_org.incoords_v + self.num_boundary_points
332+
if getattr(self.source_org, "incoords_v", None) is not None
333+
else None
334+
)
335+
incoords_w_ext = (
336+
self.source_org.incoords_w + self.num_boundary_points
337+
if getattr(self.source_org, "incoords_w", None) is not None
338+
else None
339+
)
325340
self.extended_source = fullwave.Source(
326341
p0=self.source_org.p0,
327342
coords=self.source_org.incoords + self.num_boundary_points,
328343
grid_shape=extended_grid_shape,
329344
p0_additive=self.source_org.p0_additive,
330345
coords_additive=incoords_add_ext,
346+
u0=getattr(self.source_org, "u0", None),
347+
coords_u=incoords_u_ext,
348+
v0=getattr(self.source_org, "v0", None),
349+
coords_v=incoords_v_ext,
350+
w0=getattr(self.source_org, "w0", None),
351+
coords_w=incoords_w_ext,
331352
)
332353
logger.debug("building extended source for pml...done")
333354

@@ -1434,7 +1455,7 @@ def plot(
14341455
class PMLBuilderExponentialAttenuation(PMLBuilder):
14351456
"""A class to set up PML for exponential attenuation media."""
14361457

1437-
def __init__(
1458+
def __init__( # noqa: PLR0915
14381459
self,
14391460
grid: fullwave.Grid,
14401461
medium: fullwave.Medium,
@@ -1594,12 +1615,33 @@ def __init__(
15941615
if getattr(self.source_org, "incoords_add", None) is not None
15951616
else None
15961617
)
1618+
incoords_u_ext = (
1619+
self.source_org.incoords_u + self.num_boundary_points
1620+
if getattr(self.source_org, "incoords_u", None) is not None
1621+
else None
1622+
)
1623+
incoords_v_ext = (
1624+
self.source_org.incoords_v + self.num_boundary_points
1625+
if getattr(self.source_org, "incoords_v", None) is not None
1626+
else None
1627+
)
1628+
incoords_w_ext = (
1629+
self.source_org.incoords_w + self.num_boundary_points
1630+
if getattr(self.source_org, "incoords_w", None) is not None
1631+
else None
1632+
)
15971633
self.extended_source = fullwave.Source(
15981634
p0=self.source_org.p0,
15991635
coords=self.source_org.incoords + self.num_boundary_points,
16001636
grid_shape=extended_grid_shape,
16011637
p0_additive=self.source_org.p0_additive,
16021638
coords_additive=incoords_add_ext,
1639+
u0=getattr(self.source_org, "u0", None),
1640+
coords_u=incoords_u_ext,
1641+
v0=getattr(self.source_org, "v0", None),
1642+
coords_v=incoords_v_ext,
1643+
w0=getattr(self.source_org, "w0", None),
1644+
coords_w=incoords_w_ext,
16031645
)
16041646
extended_sensor_grid_shape = tuple(
16051647
s + 2 * self.num_boundary_points for s in self.sensor_org.grid_shape

0 commit comments

Comments
 (0)