Generalized source

[joao-bapdm]

It would be helpful to have a "generalized source" that instead of acting at a single coordinate, is a field depending on the coordinate position. so instead of s = f(t), the source would be s = f(t) * g(x).

I've tried creating one source for each grid point, but the code get's very slow.

The application would be to test the convergence for the variable density case.

[jaimesouza]

Is it better than the initial conditions grid?

[joao-bapdm]

Yes.

[krober10nd]

what you're describing is already implemented

s = f(f) * g(x) 

f(t) = Ricker wavelet g(x) = Kaiser window

Anyway, you can make the the half width very large and set the call-back for f(t) return a field that's spatially variable and the size of the entire grid.

[jaimesouza]

Ok! I'd like to discuss this better later. First I gonna fix a bug in the source/receiver parallelization. I was about to change a lot in the code to fix it, but there is a simpler solution.

[joao-bapdm]

I've created a simple source class which has a different call-back function and made a horrible hack to pass the correct spacial values.

I just pass src_points and src_values directly to forward. Currently I am getting a segmentation error:

Shared object already compiled in: /home/joao/Repositories/simwave/tmp/fa343ca7ca9477abb57ce4d213b7d1a46b3d1b3e.so
double free or corruption (out)
Abortado (imagem do núcleo gravada)

I think these simple modifications should work though. If anyone wants to take a peek, I think the best way would be to add my repo as a remote and checkout this branch.

[jaimesouza]

I am gonna take a look.

[joao-bapdm]

Did I do something wrong here, @jaimesouza? I've only changed a single RickerWavelet to a MultiWavelet and passed it as an argument to wavelet in Solver. The added object is defined at line 71.

from simwave import (
    SpaceModel, TimeModel, MultiWavelet, 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 -ffast-math -Wall -std=c99 -shared'
)

# Velocity model
vel = np.zeros(shape=(512, 512), dtype=np.float32)
vel[:] = 1500.0
vel[100:] = 2000.0

# Density model
density = np.zeros(shape=(512, 512), dtype=np.float32)
density[:] = 10
density[-100:] = 20

# create the space model
space_model = SpaceModel(
    bounding_box=(0, 5120, 0, 5120),
    grid_spacing=(10, 10),
    velocity_model=vel,
    density_model=density,
    space_order=2,
    dtype=np.float32
)

# config boundary conditions
# (none,  null_dirichlet or null_neumann)
space_model.config_boundary(
    damping_length=0,
    boundary_condition=(
        "null_dirichlet", "null_dirichlet",
        "null_dirichlet", "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=[(2000, 3000), (3000, 2000)],
    window_radius=4
)

# crete the set of receivers
receiver = Receiver(
    space_model=space_model,
    coordinates=[(2560, i) for i in range(0, 5120, 10)],
    window_radius=4
)

# create a ricker wavelet with 10hz of peak frequency
ricker = RickerWavelet(10.0, time_model)
# create a multi wavelet with peak frequencies of 5Hz and 10Hz
ricker_1 = RickerWavelet(15.0, time_model)
ricker_2 = RickerWavelet(5.0, time_model)
wavelets = MultiWavelet(
    np.vstack((ricker_1.values, ricker_2.values)).T,
    time_model
)

# create the solver
solver = Solver(
    space_model=space_model,
    time_model=time_model,
    sources=source,
    receivers=receiver,
    wavelet=wavelets,
    # 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_velocity_model(space_model.density_model, file_name="density_model")
plot_wavefield(u_full[-1])
plot_shotrecord(recv)

This is the error message:

ArgumentError                             Traceback (most recent call last)
~/Repositories/simwave/mms_2D_density_v3.py in <module>
     89 
     90 # run the forward
---> 91 u_full, recv = solver.forward()
     92 
     93 print("u_full shape:", u_full.shape)

~/Repositories/simwave/simwave/kernel/frontend/solver.py in forward(self)
    127             self.receivers.interpolated_points_and_values
    128 
--> 129         u_full, recv = self._middleware.exec(
    130             operator='forward',
    131             u_full=self.u_full,

~/Repositories/simwave/simwave/kernel/backend/middleware.py in exec(self, operator, **kwargs)
    103         # run the forward operator
    104         if operator == 'forward':
--> 105             return self._exec_forward(**kwargs)
    106 
    107     def _exec_forward(self, **kwargs):

~/Repositories/simwave/simwave/kernel/backend/middleware.py in _exec_forward(self, **kwargs)
    154 
    155         # run the C forward function
--> 156         exec_time = forward(*args)
    157 
    158         print('Run forward in %f seconds.' % exec_time)

ArgumentError: argument 5: <class 'TypeError'>: array must have flags ['C_CONTIGUOUS']

[jaimesouza]

@joao-bapdm , The array is in F order apparently.

You need to convert the array to a contiguous array in memory (C order), like below:

wavelets = MultiWavelet(
    np.ascontiguousarray( np.vstack((ricker_1.values, ricker_2.values)).T ),
    time_model
)

I am gonna do this conversion in the class as default to avoid this issue.

Thanks!