Version:

Acoustic simulations in isotropic and anisotropic media

This micro-tutorial compares different parameterizations of acoustic media and demonstrates how to run VTI acoustic simulations with SalvusProject. The concept is shown for a simple 2-D layered medium, but it works just the same in 3-D by just adding an extra dimension.
Copy
import numpy as np
import os

import salvus.namespace as sn

# Run simulations on this site.
SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")

Create the Project

# Either load an existing project or create a new one.
p = sn.Project.from_domain(
    path="project",
    # Should match the volume model.
    domain=sn.domain.dim2.BoxDomain(x0=0, x1=3000, y0=-1000, y1=0),
    load_if_exists=True,
)
We start with an isotropic acoustic model with 3 layers using a custom bm model file.
%%writefile layered_model.bm
NAME         radial_model
ACOUSTIC     true
UNITS        m
COLUMNS      depth rho vp
    0.0   1000.0   1500.0
  200.0   1000.0   1500.0
  200.0   1300.0   1900.0
  500.0   1300.0   1900.0
  500.0   2200.0   3200.0
 4000.0   2200.0   3200.0
Writing layered_model.bm
Similarly, we introduce a VTI layered model, where the second and third layer have a higher wave speed in horizontal direction.
%%writefile layered_model_ani.bm
NAME         radial_model
ACOUSTIC     true
ANISOTROPIC  true
UNITS        m
COLUMNS      depth rho vpv vph
    0.0   1000.0   1500.0   1500.0
  200.0   1000.0   1500.0   1500.0
  200.0   1300.0   1900.0   2100.0
  500.0   1300.0   1900.0   2100.0
  500.0   2200.0   3200.0   3500.0
 4000.0   2200.0   3200.0   3500.0
Writing layered_model_ani.bm
The anisotropic parameters for vertical and horizontal P-wave velocities cam also represent an isotropic medium equivalent to the first example when using the same value for VPV and VPH.
%%writefile layered_model_fake_ani.bm
NAME         radial_model
ACOUSTIC     true
ANISOTROPIC  true
UNITS        m
COLUMNS      depth rho vpv vph
    0.0   1000.0   1500.0   1500.0
  200.0   1000.0   1500.0   1500.0
  200.0   1300.0   1900.0   1900.0
  500.0   1300.0   1900.0   1900.0
  500.0   2200.0   3200.0   3200.0
 4000.0   2200.0   3200.0   3200.0
Writing layered_model_fake_ani.bm
Now, we create simulation configurations for all three 1D models. Obviously, the fake anisotropic model will result in the same waveforms as the isotropic model.
for name, bm in zip(
    ["isotropic", "anisotropic", "fake-anisotropic"],
    ["layered_model.bm", "layered_model_ani.bm", "layered_model_fake_ani.bm"],
):
    p.add_to_project(
        sn.SimulationConfiguration(
            name=name,
            max_frequency_in_hertz=20.0,
            elements_per_wavelength=1.5,
            model_configuration=sn.ModelConfiguration(
                background_model=sn.model.background.one_dimensional.FromBm(
                    filename=bm, reference_datum=0.0
                ),
            ),
            event_configuration=sn.EventConfiguration(
                wavelet=sn.simple_config.stf.Ricker(center_frequency=10.0),
                waveform_simulation_configuration=sn.WaveformSimulationConfiguration(
                    end_time_in_seconds=1.5,
                ),
            ),
            absorbing_boundaries=sn.AbsorbingBoundaryParameters(
                reference_velocity=2200.0,
                reference_frequency=10.0,
                number_of_wavelengths=7,
            ),
        ),
    )


p.viz.nb.simulation_setup("anisotropic")
[2024-07-05 10:07:41,612] INFO: Creating mesh. Hang on.
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7f66e456a490>
Finally, we highlight two more model parameterizations for anisotropic acoustic material using the Thomsen parameters and a second-order tensor, which can be seen as a reduced form of the 4th-order elastic tensor.
We add two more meshes representing the equivalent isotropic and anisotropic media in this parameterization as UnstructuredMeshSimulationConfiguration.
def convert_thomsen_params_to_vti_tensor(
    dim: int,
    rho: np.ndarray,
    vp: np.ndarray,
    delta: np.ndarray,
    epsilon: np.ndarray,
):
    if dim == 2:
        model = {
            "M0": 1 / (vp**2 * rho),
            "D11": (1 + 2 * epsilon) / rho,
            "D12": (delta - epsilon) / rho,
            "D22": 1 / rho,
        }
    elif dim == 3:
        model = {
            "M0": 1 / (vp**2 * rho),
            "D11": (1 + 2 * epsilon) / rho,
            "D12": np.zeros_like(rho),
            "D13": (delta - epsilon) / rho,
            "D22": (1 + 2 * epsilon) / rho,
            "D23": (delta - epsilon) / rho,
            "D33": 1 / rho,
        }

    return model


def interpolate_vti_model_onto_mesh(mesh: sn.UnstructuredMesh):
    if "VPV" in mesh.elemental_fields.keys():
        vti_model = convert_thomsen_params_to_vti_tensor(
            dim=2,
            rho=mesh.elemental_fields["RHO"],
            vp=mesh.elemental_fields["VPV"],
            delta=0.5
            * (
                mesh.elemental_fields["VPH"] ** 2
                / mesh.elemental_fields["VPV"] ** 2
                - 1.0
            ),
            epsilon=0.5
            * (
                mesh.elemental_fields["VPH"] ** 2
                / mesh.elemental_fields["VPV"] ** 2
                - 1.0
            ),
        )

        for k, v in vti_model.items():
            mesh.elemental_fields[k] = v

        del mesh.elemental_fields["VPV"]
        del mesh.elemental_fields["VPH"]
        del mesh.elemental_fields["RHO"]

    else:
        vti_model = convert_thomsen_params_to_vti_tensor(
            dim=2,
            rho=mesh.elemental_fields["RHO"],
            vp=mesh.elemental_fields["VP"],
            delta=np.zeros_like(mesh.elemental_fields["VP"]),
            epsilon=np.zeros_like(mesh.elemental_fields["VP"]),
        )

        for k, v in vti_model.items():
            mesh.elemental_fields[k] = v

        del mesh.elemental_fields["VP"]
        del mesh.elemental_fields["RHO"]

    return mesh


mesh = p.simulations.get_mesh("isotropic")
mesh = interpolate_vti_model_onto_mesh(mesh)
p.add_to_project(
    sn.UnstructuredMeshSimulationConfiguration(
        unstructured_mesh=mesh,
        name="thomsen-isotropic",
        event_configuration=p.entities.get(
            "simulation_configuration", "anisotropic"
        ).event_configuration,
    )
)

mesh = p.simulations.get_mesh("anisotropic")
mesh = interpolate_vti_model_onto_mesh(mesh)
p.add_to_project(
    sn.UnstructuredMeshSimulationConfiguration(
        unstructured_mesh=mesh,
        name="thomsen-anisotropic",
        event_configuration=p.entities.get(
            "simulation_configuration", "anisotropic"
        ).event_configuration,
    )
)
[2024-07-05 10:07:41,912] INFO: Creating mesh. Hang on.
At the surface of the domain, we consider a single point source at the center of the horizontal extent and model a linear array of receivers with a spacing of 20 m.
src = sn.simple_config.source.cartesian.ScalarPoint2D(x=1500, y=0, f=1.0)

recs = sn.simple_config.receiver.cartesian.collections.ArrayPoint2D(
    x=np.linspace(100, 2900, 141), y=np.zeros(141), fields=["phi"]
)

p += sn.Event(event_name="event", sources=src, receivers=recs)

p.viz.nb.domain()
It is time to actually run the simulations for all 5 scenarios.
for name in [
    "isotropic",
    "anisotropic",
    "fake-anisotropic",
    "thomsen-anisotropic",
    "thomsen-isotropic",
]:
    p.simulations.launch(
        simulation_configuration=name,
        events=p.events.list(),
        site_name=SALVUS_FLOW_SITE_NAME,
        ranks_per_job=1,
    )
    p.simulations.query(block=True)
[2024-07-05 10:07:42,186] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted job_2407051007296108_df597fd5ec@local
[2024-07-05 10:07:45,067] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted job_2407051007078776_46ed3d33c5@local
[2024-07-05 10:07:47,915] INFO: Creating mesh. Hang on.
[2024-07-05 10:07:47,993] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted job_2407051007998911_ac8a3fc43d@local
[2024-07-05 10:07:50,469] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted job_2407051007477440_0c3957080c@local
[2024-07-05 10:07:53,134] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted job_2407051007144099_dafda37ae8@local
Of course, regardless of the parameterization of the medium, the waveforms of all isotropic simulations coincide. The same is true for both anisotropic parameterizations. Note, however, the difference between isotropic and anisotropic material.
p.viz.nb.waveforms(
    data=["isotropic", "fake-anisotropic", "thomsen-isotropic"],
    receiver_field="phi",
)
p.viz.nb.waveforms(
    data=["anisotropic", "thomsen-anisotropic"],
    receiver_field="phi",
)
p.viz.nb.waveforms(
    data=["isotropic", "anisotropic"],
    receiver_field="phi",
)
p.waveforms.get(data_name="isotropic", events=p.events.list())[0].plot(
    component="A",
    receiver_field="phi",
    plot_types=["shotgather"],
)
p.waveforms.get(data_name="anisotropic", events=p.events.list())[0].plot(
    component="A",
    receiver_field="phi",
    plot_types=["shotgather"],
)
PAGE CONTENTS