Mondaic - Full Waveform Solutions

So far, we've discussed different symmetry classes of elastic and acoustic materials. These symmetry classes are defined by specific mathematical transformations (e.g. rotations, reflections) that leave the material properties unchanged. For example, a hexagonal (transversely isotropic) material doesn't change its properties when rotated around the axis of symmetry.

In Salvus material, by default, the axes of the symmetry classes are aligned with the coordinate axes. However, in many cases, an actual material may be oriented arbitrarily in space, or even spatially heterogeneously oriented.:

Geological formations - Layered sedimentary rock formations frequently have tilted or folded layers
Composite materials - Engineered materials like fiber composites have specific fiber orientations
Crystalline structures - Natural crystals and minerals have preferred orientations
Fluid-filled fractures - Aligned fluid-filled cracks or fractures create directional anisotropy

Being able to specify material orientations allows modelers to represent these real-world scenarios more accurately, leading to more precise simulation results.

The material module in Salvus allows for the orientation of materials to be specified. One very important note however is that if a material like hexagonal is oriented such that the symmetry axis or planes do not anymore align with the coordinate axes, the material is no longer hexagonal in the global coordinate system. This means that after applying an orientation, materials in general will be triclinic (elastic) or orthotropic (acoustic).

For example, if we take a VTI (Vertically Transverse Isotropic) material with its symmetry axis aligned with the z-axis, and tilt it by some angle, the resulting material will need to be represented in a lower symmetry class (triclinic) when expressed in the global coordinate system.

A big caveat of using oriented materials is that the module does not easily support orienting materials for 2 dimensions. If this is a feature you are interested in, please reach out to support to get you underway on a per-case basis, and to indicate there is interest in this feature being part of the module.

from salvus.material import elastic, orientation
import numpy as np
from salvus.material._details.material import to_solver_form

np.set_printoptions(precision=2)

Salvus provides multiple ways to specify the orientation of a material. These different methods allow you to describe the orientation in terms that make the most sense for your particular application or dataset.

The AzimuthDip orientation is one of the more common ways to specify material orientation. It uses:

Azimuth: The angle (in degrees) measured from the first global axis (the x-axis) in the horizontal plane.
Dip: The angle (in degrees) measured downward from the horizontal plane.

This parameterization is particularly useful in geophysics and structural geology where field measurements are often reported in these terms.

orientation_1 = orientation.AzimuthDip.from_params(azimuth=0, dip=22.5)

orientation_1

<IPython.core.display.HTML object>

salvus.material.orientation.AzimuthDip

All orientations know how to generate a rotation matrix that can be used to rotate a material tensor (which is done internally in Salvus later):

orientation_1.generate_rotation_matrix()

array([[[ 0.92,  0.  ,  0.38],
        [ 0.  ,  1.  ,  0.  ],
        [-0.38,  0.  ,  0.92]]])

The visualization shows how the local coordinate system is oriented relative to the global coordinate system. The local z-axis (blue) is tilted at 22.5° from the global z-axis, towards the global x-axis (red).

_ = orientation.visualize_local_bases(orientation_1)

We can also specify spatially varying orientations by providing arrays for azimuth and dip:

orientation_2 = orientation.AzimuthDip.from_params(
    azimuth=0, dip=np.array([0, 11.25, 22.5, 33.75, 45.0])
)

# The semi-colon at the end of the line suppresses the output in Jupyter,
# such that the figure is not printed twice.
_ = orientation.visualize_local_bases(orientation_2)

That was different dip angles in the global x-z plane. Specifying varying azimuths would look like this:

orientation_3 = orientation.AzimuthDip.from_params(
    azimuth=np.array([0, 30, 60, 90]), dip=33
)

_ = orientation.visualize_local_bases(orientation_3)

Another way to specify orientation is using the AxesDips class, which allows you to specify:

dip_x: The dip angle of the x-axis from the horizontal plane
dip_y: The dip angle of the y-axis from the horizontal plane

This parameterization can be more intuitive when working with materials where you know the orientation of specific planes or axes relative to the global coordinate system.

orientation_4 = orientation.AxesDips.from_params(dip_x=22.5, dip_y=12.5)

orientation_4

<IPython.core.display.HTML object>

salvus.material.orientation.AxesDips

What's important to note is that the dip_x and dip_y angles don't specify the total dip of the local horizontal plane, except when one of them is 0:

orientation_5 = orientation.AxesDips.from_params(dip_x=22.5, dip_y=0)
np.allclose(
    orientation_1.generate_rotation_matrix(),
    orientation_5.generate_rotation_matrix(),
)

True

# But when dip_x and dip_y are both non-zero:
np.allclose(
    orientation_1.generate_rotation_matrix(),
    orientation_4.generate_rotation_matrix(),
)

False

This effect can be very deceptive. For relatively steep dipping planes, the visual difference is small:

orientation_5a = orientation.AxesDips.from_params(dip_x=33, dip_y=33)
orientation_5a

<IPython.core.display.HTML object>

salvus.material.orientation.AxesDips

orientation_5b = orientation.AzimuthDip.from_params(azimuth=45, dip=33)
orientation_5b

<IPython.core.display.HTML object>

salvus.material.orientation.AzimuthDip

orientation.visualize_local_bases(orientation_5a)
orientation.visualize_local_bases(orientation_5b)

<Figure size 1500x500 with 3 Axes>

This unapparent difference is an effect of the projection on the coordinate planes. When directly comparing the two orientations' rotation matrices, we see that they are indeed different:

np.allclose(
    orientation_5a.generate_rotation_matrix(),
    orientation_5b.generate_rotation_matrix(),
)

False

A final test is to convert one of the orientations into the other's system:

orientation.AzimuthDip.from_material(orientation_5a).DIP.p

42.564398528805846

orientation_5b.DIP.p

33.0

The AxesDips and AzimuthDip orientations are unfortunately not always sufficient to describe the orientation of a material. For example, if we have a material that is not hexagonal, but orthotropic, monoclinic or triclinic, we need to also keep track of the orientation of the material's horizontal rotation about the local z-axis.

Both the aforementioned orientations are not able to do this, and will simply ambiguously transform local_x and local_y to lie in the dipping plane.

A further way to specify orientation is using the AzimuthDipRake orientation, which has the same definition as AzimuthDip, but adds a rake angle. The rake angle is defined as an angle measured from the horizontal on the dipping plane (strike) for local y, or as the angle between the down direction on the dipping plane and the local x-axis. This way, when "facing" in the down direction on the dipping plane, local x is pointing in your view direction, and local y is pointing to your left. The rake angle is measured in the plane of the dip.

This is a common way to specify the orientation of faults or fractures in geophysics, where the rake angle describes the direction of slip or movement along the fault plane.

orientation_6 = orientation.AzimuthDipRake.from_params(
    azimuth=0,  # Dip towards the x-axis
    dip=22.5,
    # Local x as to make these angles with the local down, in-plane.
    rake=np.linspace(0, 30, 5),
)

orientation.visualize_local_bases(orientation_6)

<Figure size 1500x500 with 3 Axes>

Finally, we can also specify the orientation using a DirectOrthonormalBasis class. This class allows you to specify the local coordinate system using three orthonormal basis vectors. This is a more general way to specify orientation and can be useful when you have a specific coordinate system in mind or when the other parameterizations are not sufficient.

orientation_7 = orientation.DirectOrthonormalBasis.from_basis_vectors(
    local_x=np.array([0, 1.0, 0]),
    local_y=np.array([-1.0, 0.0, 0]),
    local_z=np.array([0, 0, 1.0]),
)

orientation.visualize_local_bases(orientation_7)

<Figure size 1500x500 with 3 Axes>

Note that this orientation also applies the Gram-Schmidt process to ensure that the basis vectors are orthonormal. This means that the vectors are guaranteed to be orthogonal and have unit length, but the input vectors are not necessarily respected.

orientation_8 = orientation.DirectOrthonormalBasis.from_basis_vectors(
    local_x=np.array([0, 1.0, 0]),
    local_y=np.array([-1.0, 0.0, 0]),
    local_z=np.array([0.2, 0.2, 1.0]),
    skip_orthonormalization=False,
)
orientation.visualize_local_bases(orientation_8)

orientation_8.get_basis_vectors()

(array([[ 0.  ,  0.98, -0.19]]),
 array([[-0.98,  0.04,  0.19]]),
 array([[0.2 , 0.19, 0.96]]))

Now that we've explored how to create different orientations, let's see how to apply them to materials. When we apply an orientation to a material, we're essentially rotating the material's properties to align with our specified orientation.

Here, we'll start with a VTI (Vertically Transverse Isotropic) material defined using Thomsen parameters. Thomsen parameters are commonly used in exploration geophysics to describe weak anisotropy:

epsilon (ε): Controls the difference between horizontal and vertical P-wave velocities
delta (δ): Controls the angular dependence of P-wave velocity
gamma (γ): Controls the difference between horizontal and vertical S-wave velocities

First, let's create our VTI material:

elastic_vti_material = elastic.transversely_isotropic.Thomsen.from_params(
    rho=2500, vp=5000, vs=3500, epsilon=0.23, delta=0.4, gamma=0.2
)
elastic_vti_material

<IPython.core.display.HTML object>

salvus.material.elastic.hexagonal.Thomsen

elastic_vti_material.to_tensor_components()

<IPython.core.display.HTML object>

salvus.material.elastic.hexagonal.TensorComponents

When we apply an orientation to our VTI material, we're essentially tilting its symmetry axis. In this case, we're using our AzimuthDip orientation to tilt the symmetry axis by 22.5 degrees from the vertical:

orientation_1 = orientation.AzimuthDip.from_params(azimuth=0, dip=22.5)
elastic_vti_material_oriented = elastic_vti_material.with_orientation(
    orientation_1
)

elastic_vti_material_oriented

<IPython.core.display.HTML object>

abc.salvus.material.orientation.elastic.hexagonal.Thomsen

Now let's convert our oriented material to a solver-compatible form. Salvus requires materials to be in one of its supported forms for simulation. When we apply an orientation to a material, we need to convert it to a form that the solver can understand - in this case, a triclinic tensor components material.

# to_solver_form(m=elastic_vti_material) # Won't work, nothing to do!

elastic_vti_material_oriented_sf = to_solver_form(
    m=elastic_vti_material_oriented
)

elastic_vti_material_oriented_sf

<IPython.core.display.HTML object>

salvus.material.elastic.triclinic.TensorComponents

elastic_vti_material_oriented_sf.to_wavelength_oracle()

Constant parameter
Spatially constant
Value: 3.4e+03

isinstance(
    elastic_vti_material_oriented_sf, elastic.triclinic.TensorComponents
)

True

Rotating in the symmetry plane does not break the symmetry of the tensor, but will downcast the material to triclinic. Downcasting means that the material is no longer in the original symmetry class (VTI) but is now represented in a lower symmetry class.

orientation_rotation_vti = (
    orientation.DirectOrthonormalBasis.from_basis_vectors(
        local_x=np.array([0, 1.0, 0]),
        local_y=np.array([-1.0, 0, 0]),
        local_z=np.array([0, 0, 1.0]),
    )
)

to_solver_form(elastic_vti_material.with_orientation(orientation_rotation_vti))

<IPython.core.display.HTML object>

salvus.material.elastic.triclinic.TensorComponents

Now let's create an orthotropic material using the AzimuthDipRake orientation. This material has different properties in three orthogonal directions, and we can specify its orientation using the azimuth, dip, and rake angles.

# Completely fictitious material properties
orthotropic_material = elastic.orthotropic.EngineeringConstants.from_params(
    rho=2500,
    e1=50e9,
    e2=30e9,
    e3=20e9,
    v12=0.2,
    v13=0.1,
    v23=0.3,
    g12=10e9,
    g13=5e9,
    g23=2e9,
)

orthotropic_material

<IPython.core.display.HTML object>

salvus.material.elastic.orthotropic.EngineeringConstants

# Let's visualize the entries of the stiffness tensor as well:
orthotropic_material.to_tensor_components()

<IPython.core.display.HTML object>

salvus.material.elastic.orthotropic.TensorComponents

We now first start off easy: a material with simply a 90 degree rotation about the (local) z-axis. Let's see what this does to the stiffness tensor:

orientation_orthotropic_1 = orientation.AzimuthDipRake.from_params(
    azimuth=0, dip=0, rake=90
)

orthotropic_material_oriented = orthotropic_material.with_orientation(
    orientation_orthotropic_1
)

to_solver_form(orthotropic_material_oriented)

<IPython.core.display.HTML object>

salvus.material.elastic.triclinic.TensorComponents

This resulting material has the entries for C11 and C22 swapped, C13 and C23 swapped, and C44 and C55 swapped. Additionally, some numerical noise has entered some of the off-diagonal terms. This is not an issue when simulating with Salvus, but can be manually muted if desired.

To map this material back to an orthotropic form (and mute numerical noise), Salvus will require you to forcefully discard information, as strictly speaking, the material is now triclinic.

elastic.orthotropic.EngineeringConstants.from_material(
    to_solver_form(orthotropic_material_oriented),
    reduction_method="remove-components",
)

<IPython.core.display.HTML object>

salvus.material.elastic.orthotropic.EngineeringConstants

Also note that remove-components only works if the entries in the to-be removed components are actually close to zero. This is apparent when we try the same approach with a less than 90 degree rotation, under which the material doesn't remain orthotropic in the global coordinate system:

orthotropic_material_oriented_45 = orthotropic_material.with_orientation(
    orientation.AzimuthDipRake.from_params(azimuth=0, dip=0, rake=45)
)
to_solver_form(orthotropic_material_oriented_45)

<IPython.core.display.HTML object>

salvus.material.elastic.triclinic.TensorComponents

try:
    elastic.orthotropic.EngineeringConstants.from_material(
        to_solver_form(orthotropic_material_oriented_45),
        reduction_method="remove-components",
    )
except TypeError as e:
    print(e)

Passed material of type <class 'salvus.material.elastic.triclinic.TensorComponents'> can't be treated as <class 'salvus.material.elastic.orthotropic.TensorComponents'>, as it does not respect the symmetry. The component C16 should be zero, but it is not.

The "force" method will work, but information will be lost and the material will be altered:

elastic.orthotropic.EngineeringConstants.from_material(
    to_solver_form(orthotropic_material_oriented_45), reduction_method="force"
)

<IPython.core.display.HTML object>

salvus.material.elastic.orthotropic.EngineeringConstants

Finally, let's see how an orthotropic material can map into a fully triclinic material.

orthotropic_material_oriented_dramatic = orthotropic_material.with_orientation(
    orientation.AzimuthDipRake.from_params(azimuth=22.5, dip=30, rake=45)
)
to_solver_form(orthotropic_material_oriented_dramatic)

<IPython.core.display.HTML object>

salvus.material.elastic.triclinic.TensorComponents

This even produces negative moduli in the off-diagonal terms!

In the final tutorial of this series, we will explore how to set up a simulation using the materials we've created.

Part 2 - Anisotropy and elliptic coefficients

Part 4 - Simulating with Materials and Exotic Waveforms from Anisotropy

© Mondaic AG (2025)

Site Map Contact Us

Academic License Agreement

PAGE CONTENTS