Mondaic
PricingAbout UsDocumentation
Mondaic
Version:
No results found.

This documentation is not for the latest stable Salvus version.

Meshing cylindrical domains

Several 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.
Copy
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.0

Mesh creation

The 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     1500.0       0.0
    0.1   1000.0     1500.0       0.0
    0.1   7710.0     5890.0    3240.0
    0.12  7710.0     5890.0    3240.0
    0.12  1000.0     1500.0       0.0
    0.5   1000.0     1500.0       0.0
Writing radial_model.bm

High-order shape transformations

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.unstructured_mesh.UnstructuredMesh at 0x7f59080fd9d0>
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.unstructured_mesh.UnstructuredMesh at 0x7f5827d35c10>
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.unstructured_mesh.UnstructuredMesh at 0x7f58217f5110>
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.unstructured_mesh.UnstructuredMesh at 0x7f5829dfca90>
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 0x7f58254ecfd0>
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-03-15 09:15:45,203] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted job_2403150915462863_2c6ed327ce@local

True
p.viz.shotgather(["sim"], "event", "phi", "A")
<Axes: title={'center': 'Event: event | field: phi | component: A'}, xlabel='Receiver Index', ylabel='Time [s]'>
Previous
Mesh-to-mesh interpolation
Next
Data-adaptive mesh masking
We use cookies to analyze our traffic, but won't track you until explicit consent is given. Using our contact form will create relevant cookies and transmit data to third party services to keep track of your request. Read more about how we use personal information in our Privacy Policy.
PAGE CONTENTS