# Simwave

`Simwave`

is a Python package to simulate the propagation of the constant or variable density acoustic wave in an isotropic 2D/3D medium using the finite difference method. Finite difference kernels of aribtrary spatial order (up to 20th order) are written in C for performance and compiled at run time. These kernels are called via a user-friendly Python interface for easy integration with several scientific and engineering libraries to, for example perform full-waveform inversion.

## Installation:

For installation, `simwave`

needs only scipy, numpy, and segyio. See `requirements.txt`

. If you wish to plot, then `matplotlib`

is additionally required. `simwave`

compiles finite difference stencils at run time in C for performance and thus requires a working C compiler.

`git clone https://github.com/HPCSys-Lab/simwave.git`

`cd simwave`

`pip3 install -e .`

## Contributing

All contributions are welcome.

To contribute to the software:

- Fork the repository.
- Clone the forked repository, add your contributions and push the changes to your fork.
- Create a Pull request

Before creating the pull request, make sure that the tests pass by running

```
pytest
```

Some things that will increase the chance that your pull request is accepted:

- Write tests.
- Add Python docstrings that follow the Sphinx.
- Write good commit and pull request messages.

# Problems?

If something isn't working as it should or you'd like to recommend a new addition/feature to the software, please let us know by starting an issue through the issues tab. I'll try to get to it as soon as possible.

# Examples

Simulation with `simwave`

is simple and can be accomplished in a dozen or so lines of Python! Jupyter notebooks with tutorials can be found here here.

Here we show how to simulate the constant density acoustic wave equation on a simple two layer velocity model.

```
from simwave import (
SpaceModel, TimeModel, RickerWavelet, Solver, Compiler,
Receiver, Source, plot_wavefield, plot_shotrecord, plot_velocity_model
)
import numpy as np
# set compiler options
# available language options: c (sequential) or cpu_openmp (parallel CPU)
compiler = Compiler(
cc='gcc',
language='cpu_openmp',
cflags='-O3 -fPIC -Wall -std=c99 -shared'
)
# Velocity model
vel = np.zeros(shape=(512, 512), dtype=np.float32)
vel[:] = 1500.0
vel[100:] = 2000.0
# create the space model
space_model = SpaceModel(
bounding_box=(0, 5120, 0, 5120),
grid_spacing=(10, 10),
velocity_model=vel,
space_order=4,
dtype=np.float32
)
# config boundary conditions
# (none, null_dirichlet or null_neumann)
space_model.config_boundary(
damping_length=0,
boundary_condition=(
"null_neumann", "null_dirichlet",
"none", "null_dirichlet"
),
damping_polynomial_degree=3,
damping_alpha=0.001
)
# create the time model
time_model = TimeModel(
space_model=space_model,
tf=1.0,
saving_stride=0
)
# create the set of sources
source = Source(
space_model,
coordinates=[(2560, 2560)],
window_radius=1
)
# crete the set of receivers
receiver = Receiver(
space_model=space_model,
coordinates=[(2560, i) for i in range(0, 5112, 10)],
window_radius=1
)
# create a ricker wavelet with 10hz of peak frequency
ricker = RickerWavelet(10.0, time_model)
# create the solver
solver = Solver(
space_model=space_model,
time_model=time_model,
sources=source,
receivers=receiver,
wavelet=ricker,
compiler=compiler
)
# run the forward
u_full, recv = solver.forward()
print("u_full shape:", u_full.shape)
plot_velocity_model(space_model.velocity_model)
plot_wavefield(u_full[-1])
plot_shotrecord(recv)
```