%matplotlib inline
# This notebook will use this variable to determine which
# remote site to run on.
import os
SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
import matplotlib.pyplot as plt
import numpy as np
import pathlib
import time
import xarray as xr
import salvus.namespace as sn
20 x 20 cm
centered around the origin and insert an object with two spherical inclusions that mimics a breast phantom in a USCT acquisition. The model for density (RHO
) and speed of sound (VP
) is created from a structured grid.def get_spherical_inclusion():
nx, ny = 200, 200
x = np.linspace(-0.1, +0.1, nx)
y = np.linspace(-0.1, +0.1, nx)
xx, yy = np.meshgrid(x, y, indexing="ij")
# Add 3 spherical inclusions
vp = 1500.0 * np.ones_like(xx)
rho = 990.0 * np.ones_like(xx)
mask = np.sqrt(xx**2 + yy**2) < 0.05
vp[mask] = 1480.0
rho[mask] = 975.0
mask = np.sqrt(xx**2 + (yy - 0.025) ** 2) < 0.015
vp[mask] = 1555.0
rho[mask] = 1040.0
mask = np.sqrt(xx**2 + (yy + 0.025) ** 2) < 0.015
vp[mask] = 1460.0
rho[mask] = 1010.0
ds = xr.Dataset(
data_vars={
"vp": (["x", "y"], vp),
"rho": (["x", "y"], rho),
},
coords={"x": x, "y": y},
)
return ds
true_model = get_spherical_inclusion()
# Plot the xarray dataset.
plt.figure(figsize=(16, 6))
plt.subplot(121)
true_model.vp.T.plot()
plt.subplot(122)
true_model.rho.T.plot()
plt.suptitle("Model with spherical inclusions")
plt.show()
true_model
.# Uncomment the following line to delete a
# potentially existing project for a fresh start
# !rm -rf project
vm = sn.model.volume.cartesian.GenericModel(name="true_model", data=true_model)
p = sn.Project.from_volume_model(
path="project",
volume_model=vm,
load_if_exists=True,
)
simple_config
has a few built-in options to create lists of sources and receivers, which we want to make use of. By defining the two rings below - one for sources and one for the receivers, we can readily create an EventCollection
, i.e., a number of experiments characterized by the locations of the emitting and the receiving transducers.srcs = sn.simple_config.source.cartesian.collections.ScalarPoint2DRing(
x=0, y=0, radius=0.09, count=5, f=1.0
)
for _i, s in enumerate(srcs._sources):
all_recs = sn.simple_config.receiver.cartesian.collections.RingPoint2D(
x=0, y=0, radius=0.09, count=100, fields=["phi"]
)
recs = [
r
for r in all_recs._receivers
if np.sqrt(
(s.location[0] - r.location[0]) ** 2
+ (s.location[1] - r.location[1]) ** 2
)
> 0.03
]
p += sn.EventCollection.from_sources(
sources=[s], receivers=recs, event_name_starting_index=_i
)
1.5 ms
(and later on also the same source time function) for each of them. Here, we also fix the time step to 400 ns
. While this is not strictly necessary and could be automatically detected, of course, it will simplify the data comparison across different simulations.wsc = sn.WaveformSimulationConfiguration(end_time_in_seconds=0.00015)
wsc.physics.wave_equation.time_step_in_seconds = 4.0e-7
50 kHz
and a mesh accurate up to 100 kHz
.ec = sn.EventConfiguration(
waveform_simulation_configuration=wsc,
wavelet=sn.simple_config.stf.Ricker(center_frequency=50000.0),
)
p += sn.SimulationConfiguration(
name="true_model_100kHz",
#
# Settings that influence the mesh.
elements_per_wavelength=2,
tensor_order=4,
max_frequency_in_hertz=100000.0,
#
model_configuration=sn.ModelConfiguration(
background_model=None, volume_models="true_model"
),
# Potentially event dependent settings.
event_configuration=ec,
)
p.simulations.launch(
simulation_configuration="true_model_100kHz",
events=p.events.list(),
site_name=SALVUS_FLOW_SITE_NAME,
ranks_per_job=1,
)
[2024-11-06 22:30:50,729] INFO: Creating mesh. Hang on. [2024-11-06 22:30:50,799] INFO: Submitting job array with 5 jobs ...
5
p.simulations.query(block=True)
True
p.simulations.get_mesh("true_model_100kHz")