import salvus.namespace as sn
import numpy as np
import os
SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
PROJECT_DIR = "project"
v_min = 1500.0
elements_per_wavelength = 4.0
frequency_in_hertz = 5e3
length_in_meters = 1.0
%%writefile radial_model.bm
NAME radial_model
UNITS m
COLUMNS radius rho vp vs
0.0 1000.0 2000.0 0.0
0.1 1000.0 2000.0 0.0
0.1 7710.0 5890.0 3240.0
0.12 7710.0 5890.0 3240.0
0.12 1000.0 2000.0 0.0
0.5 1000.0 2000.0 0.0
Writing radial_model.bm
tensor_order
within the Salvus ecosystem, and applies to both the polynomial degree of the model as well as the shape transformation.tensor_order
of 1 and 4, respectively.mc = sn.simple_mesh.Circular2D()
mc.basic.model = "radial_model.bm"
mc.basic.min_period_in_seconds = 1 / 2e3
mc.basic.elements_per_wavelength = elements_per_wavelength
mc.validate()
mc.create_mesh()
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7f0c5053d6d0>
mc.advanced.tensor_order = 4
mesh = mc.create_mesh()
mesh
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7f0d1e70b8d0>
tensor_order
and only depend on the polynomial degree of the spectral-element shape functions (which default to 4). Hence, by changing the tensor order to 4, we obtain a much more accurate simulation at the same cost.tensor_order
of 1 as input, but we change the tensor order later on, to ensure a smooth interpolation of the curved shape.mc = sn.simple_mesh.Circular2D()
mc.basic.model = "radial_model.bm"
mc.basic.min_period_in_seconds = 1 / frequency_in_hertz
mc.basic.elements_per_wavelength = elements_per_wavelength
mc.validate()
mesh_2d = mc.create_mesh()
mesh_2d
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7f0c4a2b6e90>
r1
in radial direction, z0
at the bottom, and z1
at the top.def edge_length(
v_min: float, elements_per_wavelength: float, frequency_in_hertz: float
):
return v_min / (frequency_in_hertz * elements_per_wavelength)
offsets = np.linspace(
0.0,
length_in_meters,
int(
np.ceil(
length_in_meters
/ edge_length(v_min, elements_per_wavelength, frequency_in_hertz)
)
),
endpoint=True,
)
mesh = mesh_2d.extrude(offsets)
mesh.find_side_sets("cylindrical")
# Change tensor order to get a proper circular shape.
mesh.change_tensor_order(4, interpolation_mode="cylindrical")
radius_2d = mesh.get_element_centroid()[:, :2]
mask1 = np.linalg.norm(radius_2d, axis=1) > 0.1
mask2 = np.linalg.norm(radius_2d, axis=1) < 0.12
mask = mask1 & mask2
vp_water = 1500.0
rho_water = 1000.0
vp_steel = 5890.0
vs_steel = 3240.0
rho_steel = 7710.0
vp = vp_water * np.ones_like(mesh.connectivity)
vp[mask] = vp_steel
mesh.attach_field("VP", vp)
vs = np.zeros_like(mesh.connectivity)
vs[mask] = vs_steel
mesh.attach_field("VS", vs)
rho = rho_water * np.ones_like(mesh.connectivity)
rho[mask] = rho_steel
mesh.attach_field("RHO", rho)
fluid = np.ones(mesh.nelem)
fluid[mask] = 0.0
mesh.attach_field("fluid", fluid)
mesh
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7f0c4a2801d0>
d = sn.domain.dim3.BoxDomain(x0=-0.5, x1=0.5, y0=-0.5, y1=0.5, z0=0.0, z1=1.0)
p = sn.Project.from_domain(path=PROJECT_DIR, domain=d, load_if_exists=True)
p.add_to_project(
sn.Event(
event_name="event",
sources=[
sn.simple_config.source.cartesian.ScalarPoint3D(
x=0.0,
y=0.0,
z=0.5,
f=1.0,
)
],
receivers=[
sn.simple_config.receiver.cartesian.Point3D(
x=0.0,
y=0.0,
z=z,
network_code="AA",
station_code=f"{_i:02}",
fields=["phi"],
)
for _i, z in enumerate(np.linspace(0.0, 1.0, 20))
],
)
)
UnstructuredMeshSimulationConfiguration
, which is just a special case of the SimulationConfiguration
. After adding the configuration to the project, it can be used in the same way as a configuration that uses the built-in meshing routines.if "sim" not in p.simulations.list():
p.add_to_project(
sn.UnstructuredMeshSimulationConfiguration(
name="sim",
unstructured_mesh=mesh,
event_configuration=sn.EventConfiguration(
wavelet=sn.simple_config.stf.Ricker(
center_frequency=frequency_in_hertz
),
waveform_simulation_configuration=sn.WaveformSimulationConfiguration(
end_time_in_seconds=0.00035,
),
),
),
)
p.viz.nb.simulation_setup("sim", "event")
<salvus.flow.simple_config.simulation.waveform.Waveform object at 0x7f0d34278bd0>
p.simulations.launch(
"sim",
p.events.list(),
site_name=SALVUS_FLOW_SITE_NAME,
ranks_per_job=2,
)
p.simulations.query(block=True)
[2024-07-05 10:05:32,286] INFO: Submitting job ... Uploading 1 files... 🚀 Submitted job_2407051005395007_297a7baff3@local
True
p.viz.shotgather(["sim"], "event", "phi", "A")
<Axes: title={'center': 'Event: event | field: phi | component: A'}, xlabel='Receiver Index', ylabel='Time [s]'>