# Variable used in the notebook to determine which site
# is used to run the simulations.
import os
SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
PROJECT_DIR = "project"
import time
import pathlib
import matplotlib.pyplot as plt
import numpy as np
import obspy
import xarray as xr
from salvus import namespace as sn
# Define a spherical chunk domain.
lat_c, lon_c = -27.5, 135
lat_e, lon_e = 25.5, 42.0
# Uncomment the following line to start from scratch
# !rm -rf project
d = sn.domain.dim3.SphericalChunkDomain(
lat_center=lat_c,
lat_extent=lat_e,
lon_center=lon_c,
lon_extent=lon_e,
radius_in_meter=6371e3,
)
p = sn.Project.from_domain(path=PROJECT_DIR, domain=d, load_if_exists=True)
# Add an event by parsing some prepared data.
e = sn.Event(
sources=sn.simple_config.source.seismology.parse(
"data/event.txt", dimensions=3
),
receivers=sn.simple_config.receiver.seismology.parse(
obspy.read_inventory("data/stations.txt"),
dimensions=3,
fields=["displacement"],
),
)
p += e
lat = np.linspace(lat_c - lat_e / 2.0, lat_c + lat_e / 2.0, 50)
lon = np.linspace(lon_c - lon_e / 2.0, lon_c + lon_e / 2.0, 50)
depth = np.linspace(0, 1000.0, 40)
_, lon_grid, lat_grid = np.meshgrid(depth, lon, lat, indexing="ij")
# Distance from center in degree.
d = ((lon_grid - lon_c) ** 2 + (lat_grid - lat_c) ** 2) ** 0.5
# Apply Gaussian
sigma = 5.0
d = 1.0 / (sigma * np.sqrt(2 * np.pi)) * np.exp(-0.5 * (d / sigma) ** 2)
# Normalize to 10 % deviation
d /= d.max()
d *= 10.0
ds = xr.Dataset(
data_vars={
"vs": (["depth", "longitude", "latitude"], d),
"vp": (["depth", "longitude", "latitude"], d),
},
coords={"depth": depth, "latitude": lat, "longitude": lon},
)
ds.vs.attrs["units"] = "%"
ds.vp.attrs["units"] = "%"
# Same for the coordinate axes.
ds.depth.attrs["units"] = "km"
ds.latitude.attrs["units"] = "degrees_north"
ds.longitude.attrs["units"] = "degrees_east"
ds.attrs["geospatial_lat_units"] = "degrees_north"
ds.attrs["geospatial_lon_units"] = "degrees_east"
ds.attrs["geospatial_vertical_units"] = "km"
ds.vs.isel(depth=0).T.plot()
plt.show()
ec = sn.EventConfiguration(
wavelet=sn.simple_config.stf.GaussianRate(half_duration_in_seconds=100.0),
waveform_simulation_configuration=sn.WaveformSimulationConfiguration(
end_time_in_seconds=1800.0
),
)
mc = sn.ModelConfiguration(background_model="prem_iso_no_crust")
p += sn.SimulationConfiguration(
tensor_order=1,
name="initial_model",
elements_per_wavelength=1.0,
min_period_in_seconds=200,
max_depth_in_meters=500e3,
model_configuration=sn.ModelConfiguration(
background_model="prem_iso_no_crust",
),
event_configuration=ec,
)
p += sn.model.volume.seismology.GenericModel(name="vs_blob", data=ds)
p += sn.SimulationConfiguration(
tensor_order=1,
name="target_model",
elements_per_wavelength=1.0,
min_period_in_seconds=200,
max_depth_in_meters=500e3,
model_configuration=sn.ModelConfiguration(
background_model="prem_iso_no_crust", volume_models=["vs_blob"]
),
event_configuration=ec,
)
p.viz.nb.simulation_setup("target_model", events=p.events.list())
[2024-11-06 22:38:39,192] INFO: Creating mesh. Hang on.
Interpolating model: vs_blob.