MondaicSeveral applications require a cylindrical geometry, for instance, to model boreholes, pipelines or the tissue/bone structure of the human limbs.
In the simplest case, a cylindrical medium can be described by a 1D radial model, which is then extruded in the third dimension.
SalvusProject currently has no built-in support for such meshes, however, we can utilize external meshing tools. This mini tutorial shows how to create such meshes and how to add them to a project.
import salvus.namespace as sn
import numpy as np
import os
SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
PROJECT_DIR = "project"In this example, we model a cylindrical domain that mimics a steel pipe filled and surrounded by water. We use a dominant frequency of 5 kHz and a vertical extent of 1 m.
v_min = 1500.0
elements_per_wavelength = 4.0
frequency_in_hertz = 5e3
length_in_meters = 1.0The mesh generation proceeds in two steps. First, we create a disk-shaped 2D mesh that includes the radial model. In a second step, the mesh is extruded along the longitudinal axis. The radial profile can be specified as a bm file, which contains the material model as well as interfaces.
We consider a steel pipe with an inner diameter of 20 cm and a thickness of 2 cm, which is submerged in water. The mesh generation algorithm will automatically extract fluid and solid parts of the domain, and adjust the element size to the slowest velocity.
%%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.0Writing radial_model.bm
Salvus supports different polynomial orders for the shape transformations of the quadrilateral or hexahedral elements. This is particularly useful when meshing curved domain, where a piecewise (bi-)linear transformation could result in a poorly approximated shape. This parameter is called
tensor_order within the Salvus ecosystem, and applies to both the polynomial degree of the model as well as the shape transformation.To demonstrate the difference, we will take a quick detour and generate two meshes for a frequency of 2 kHz and with
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 0x7fa92f39c790>
We notice that the coarse mesh poorly models the concentric shape of the 1D model. Note that the element-size is chosen according to resolve the desired frequency at the specified number of elements per wavelength. However, the predominant numerical error is introduced by the poor approximation of the geometry.
When changing the tensor order to 4, the edges of the elements are defined by 4th-order Lagrange polynomials, which results in a much better shape.
mc.advanced.tensor_order = 4
mesh = mc.create_mesh()
mesh
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7fa91c1b0e10>
It is important to note that the number of elements is the same in both cases. Furthermore, the computational costs of the simulation are independent of the
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.Now, we switch back to the actual domain and increase the dominant frequency to 5 kHz. Note that the second step requires a mesh with a
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 0x7fa91c26ff10>
Next, we extrude the mesh along the z-axis. The horizontal mesh is kept fixed, and we only need to determine the edge length in longitudinal direction, which can be computed as a function of the minimum velocity (sound speed in water), frequency and elements per wavelength.
After extruding the mesh, we need to detect the side sets of the cylindrical mesh, which are labeled as
r1 in radial direction, z0 at the bottom, and z1 at the top.As a final step, we change the tensor order to 4.
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")The extrusion acts only on the mesh topology, but does not add the material model. We already ensured though that the outline of the steel pipe is aligned with the mesh, and we can easily attach the material model based on the 2D radius.
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 0x7fa91c255c50>
This concludes the meshing procedure. Possible extensions could include adding a deformation to the radial profile, or considering a 3D velocity structure, but we leave this to another tutorial.
Now that the mesh has been created, we start a project to make use of the simulation and visualization tools. Because SalvusProject does not yet natively support cylindrical domains, we specify a 3-D box domain enclosing the cylindrical domain.
For simplicity, we add a single event in the center of the mesh and consider a line of receivers along the longitudinal axis.
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))
],
)
)
Externally built meshes can be used in any project using an
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 0x7fa91a399dd0>
To conclude the tutorial, we launch a simple simulation and plot a shotgather of the traces.
p.simulations.launch(
"sim",
p.events.list(),
site_name=SALVUS_FLOW_SITE_NAME,
ranks_per_job=2,
)
p.simulations.query(block=True)[2024-11-20 13:17:25,020] INFO: Submitting job ... Uploading 1 files... 🚀 Submitted job_2411201317106539_6fe34545dd@local
True
p.viz.shotgather(["sim"], "event", "phi", "A")<Axes: title={'center': 'Event: event | field: phi | component: A'}, xlabel='Receiver Index', ylabel='Time [s]'>
PAGE CONTENTS