Minor suggestions for `compiler` class

[krober10nd]

1. Even though I specified the language=c compilation fails because the code tries to import omp.h header despite it not being used. Perhaps this could be avoided by putting this include in a pragma?
2. In the case that the compilation failed, you can catch this error and write a helpful error message to the screen informing the user what compiler language's are supported.

[krober10nd]

3. The hash used to save the shared object should detect if the compilation options changed and force a recompilation option if they are different. Otherwise it seems to not recompile despite changing the options.

[jaimesouza]

@krober10nd Do you have an example where the options change, but simwave does not recompile the C code? It should recompile if any of the following parameters changes:

• language
• cc
• cflags
• dimension
• operator (current only forward)
• density (constant or variable density)
• float precision
• C file code

[krober10nd]

Yes, the examples Joao posted in our slack channel. Specifically I was running (and re-running) the loop_in_time example. I switched the -O0 flag with the -O2 flag and it didn't recompile and more concerning, sometimes re-running it several times in a sequence back-to-back produced completely different results.

[jaimesouza]

I will take a look

[jaimesouza]

Fixing 1: https://github.com/HPCSys-Lab/simwave/commit/60cadef9ea970c54d8741786ebf1a93a9c766f00

[jaimesouza]

Regarding 3, as shown in the images below, the hash (also name of the shared object) chages when I change the flags.

Perhaps, you compiled the code with the same flags before and the shared object was already created in the tmp folder, located in your working dir.

Every time a new shared object is created, it shows the compilation command.

Compilation command: ...

Otherwise, it's shown:

Shared object already compiled in: ...

Can you try your example again and confirm if the behavior is as expected, please?

[jaimesouza]

more concerning, sometimes re-running it several times in a sequence back-to-back produced completely different results.

This can be a serious problem. Do you mean, running the simulation (with the same parameters) over and over in a loop, for examplo, can lead to differente results?

[krober10nd]

Yes, I see now, but perhaps there was a hash collision :P

Anyway, yes, both João and I observed this behavior from time to time. Just keep running the acoustic tests on the slack channel.

[jaimesouza]

@krober10nd or @joao-bapdm I could not find this inconsistent behavior when running several times in a sequence. Could you post here an example (input and output) where this problem accurs, please?

[joao-bapdm]

I am able to get these plots with

import numpy as np
import matplotlib.pyplot as plt

from simwave import (
    SpaceModel, TimeModel, RickerWavelet, MultiWavelet, Wavelet, Solver, Compiler,
    Receiver, Source, plot_wavefield, plot_shotrecord, plot_velocity_model
)

def grid(Lx, Lz, h):
    nz = int((Lz) / h + 1)
    nx = int((Lx) / h + 1)
    zcoords = np.linspace(0, Lz, nz)
    xcoords = np.linspace(0, Lx, nx)
    return np.meshgrid(zcoords, xcoords)

def ansatz(Lx, Lz, h, freq, tp):
    # Params
    kz = np.pi / Lx
    kx = np.pi / Lz
    meshgrid = grid(Lx, Lz, h)
    omega = 2*np.pi*freq
    # Solution
    space_solution = np.sin(kx*meshgrid[0]) * np.sin(kz*meshgrid[1])
    time_solution = np.cos(omega*tp)
    return (space_solution[:,:,None] * time_solution).T

def initial_shape(space_model, kx, kz):
    nbl_pad_width = space_model.nbl_pad_width
    halo_pad_width = space_model.halo_pad_width
    bbox = space_model.bounding_box
    # values
    xmin, xmax = bbox[0: 2]
    zmin, zmax = bbox[2: 4]
    x_coords = np.linspace(xmin, xmax, space_model.shape[0])
    z_coords = np.linspace(zmin, zmax, space_model.shape[1])
    X, Z = np.meshgrid(x_coords, z_coords)

    return np.sin(kx * X) * np.sin(kz * Z)

class MySource(Source):

    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.count = 1

    @property
    def count(self):
        """Return the number of sources/receives."""
        return self._count

    @count.setter
    def count(self, count):
        """Return the number of sources/receives."""
        self._count = count

    @property
    def interpolated_points_and_values(self):
        """ Calculate the amplitude value at all points"""

        points = np.array([], dtype=np.uint)
        values = np.array([], dtype=self.space_model.dtype)
        offset = np.array([0], dtype=np.uint)

        nbl_pad_width = self.space_model.nbl_pad_width
        halo_pad_width = self.space_model.halo_pad_width
        dimension = self.space_model.dimension
        bbox = self.space_model.bounding_box

        # senoid window
        for dim, dim_length in enumerate(self.space_model.shape):

            # points
            lpad = halo_pad_width[dim][0] + nbl_pad_width[dim][0]
            rpad = halo_pad_width[dim][1] + nbl_pad_width[dim][1] + dim_length
            points = np.append(points, np.array([lpad, rpad], dtype=np.uint))
            # values
            coor_min, coor_max = bbox[dim * dimension:2 + dim * dimension]
            coordinates = np.linspace(coor_min, coor_max, dim_length)
            amplitudes = np.sin(np.pi / (coor_max - coor_min) * coordinates)
            values = np.append(values, amplitudes.astype(values.dtype))
        # offset
        offset = np.append(offset, np.array([values.size], dtype=np.uint))

        # add window for density term

        # points
        lpad_x = halo_pad_width[0][0] + nbl_pad_width[0][0]
        rpad_x = halo_pad_width[0][1] + nbl_pad_width[0][1] + self.space_model.shape[0]
        points = np.append(points, np.array([lpad_x, rpad_x], dtype=np.uint))

        lpad_z = halo_pad_width[1][0] + nbl_pad_width[1][0]
        rpad_z = halo_pad_width[1][1] + nbl_pad_width[1][1] + self.space_model.shape[1]
        points = np.append(points, np.array([lpad_z, rpad_z], dtype=np.uint))

        xmin, xmax = bbox[0: 2]
        zmin, zmax = bbox[2: 4]
        kx = np.pi / (xmax - xmin)
        kz = np.pi / (zmax - zmin)
        x_coordinates = np.linspace(xmin, xmax, self.space_model.shape[0])
        z_coordinates = np.linspace(zmin, zmax, self.space_model.shape[1])
        X, Z = np.meshgrid(z_coordinates, z_coordinates)

        # (1 / rho) * drho_dz * dphi_dz
        x_amplitudes = np.sin(kx * x_coordinates)
        drho_dz = np.gradient(dens[:, 256], z_coordinates).mean()
        z_amplitudes = kz * drho_dz * np.cos(kz * z_coordinates) / dens[:,256]

        values = np.append(values, x_amplitudes.astype(values.dtype))
        values = np.append(values, z_amplitudes.astype(values.dtype))

        offset = np.append(offset, np.array([values.size], dtype=np.uint))

        return points, values, offset

def cosenoid(time_model, amplitude, omega=2*np.pi*1):
    """Function that generates the signal satisfying s(0) = ds/dt(0) = 0"""
    return amplitude * np.cos(omega * time_model.time_values)

def create_models(Lx, Lz, vel, dens, tf, h, p):
    # create the space model
    space_model = SpaceModel(
        bounding_box=(0, Lx, 0, Lz),
        grid_spacing=(h, h),
        velocity_model=vel,
        density_model=dens,
        space_order=p,
        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=tf,
        saving_stride=1
    )

    return space_model, time_model

def geometry_acquisition(space_model, sources=None, receivers=None):
    if sources is None:
        sources=[(2560, 2560)]
    if receivers is None:
        receivers=[(2560, i) for i in range(0, 5120, 10)]

    # create the set of sources
    # source = Source(
    #     space_model,
    #     coordinates=sources,
    #     window_radius=4
    # )

    # personalized source
    my_source = MySource(space_model, coordinates=[])

    # crete the set of receivers
    receiver = Receiver(
        space_model=space_model,
        coordinates=receivers,
        window_radius=4
    )

    return my_source, receiver

def build_solver(space_model, time_model, vp, kx, kz, omega):

    # geometry acquisition
    source, receiver = geometry_acquisition(space_model)

    # create a cosenoid wavelet
    amplitude = 1 * ((kx**2+kz**2) - omega**2/vp**2)
    # amplitude = 0
    wavelet1 = Wavelet(cosenoid, time_model=time_model, amplitude=amplitude, omega=omega)
    wavelet2 = Wavelet(cosenoid, time_model=time_model, amplitude=1, omega=omega)
    # wavelets = MultiWavelet(np.vstack((wavelet1.values,)).T, time_model)
    wavelets = MultiWavelet(np.vstack((wavelet1.values, wavelet2.values)).T, time_model)

    source.count = wavelets.values.shape[-1]

    # set compiler options
    compiler = Compiler( cc='gcc', language='cpu_openmp', cflags='-O3 -fPIC -ffast-math -Wall -std=c99 -shared')

    # create the solver
    solver = Solver(
        space_model=space_model,
        time_model=time_model,
        sources=source,
        receivers=receiver,
        wavelet=wavelets,
        compiler=compiler
    )

    return solver

if __name__ == "__main__":

    # problem params
    Lx, Lz = 5120, 5120
    tf = 2 * 1.5
    freq = 5
    vp = 1500
    rho = 1
    rho_z = 1e6
    rho_x = 0

    # discretization params
    n = 513    
    h = 10
    p = 4

    # harmonic parameters
    kx = np.pi / Lx
    kz = np.pi / Lz
    omega = 2*np.pi*freq

    # velocity model
    vel = vp * np.ones(shape=(n, n), dtype=np.float32)

    # density model
    dens = rho * np.ones(shape=(n, n), dtype=np.float32)
    for depth, dens_values in enumerate(dens):
        dens_values[:] = rho + rho_z * (dens.shape[0] - depth) / dens.shape[0] 

    # discretization
    space_model, time_model = create_models(Lx, Lz, vel, dens, tf, h, p)

    # initial grid
    initial_grid = initial_shape(space_model, kx, kz)

    # get solver
    solver = build_solver(space_model, time_model, vp, kx, kz, omega)

    # run the forward
    u_full, recv = solver.forward(
        [initial_grid * np.cos(-1 * time_model.dt * omega),
         initial_grid * np.cos(-0 * time_model.dt * omega)
        ]
    )

    # plot stuff
    print("u_full shape:", u_full.shape)
    plot_wavefield(u_full[-1])
    plot_shotrecord(recv)

    # Get analytical solution
    u_analytic = ansatz(Lx, Lz, h, freq, time_model.time_values)
    # plot comparison
    numerical_values = u_full[:, 256, 256]
    analytic_values = u_analytic[:, 256, 256]
    time_values = time_model.time_values
    plt.plot(time_values, analytic_values, time_values, numerical_values)
    plt.legend(["analytical", "numerical"])
    plt.title("Comparison between analytical and numerical result (simwave)")
    plt.xlabel("time")
    plt.show()

[jaimesouza]

@joao-bapdm, the source interpolated_points are all inclusive.

So [ 2 515 2 515 2 515 2 515] should be [ 2 514 2 514 2 514 2 514].

The problem occurs because the application is accessing an out of bound memory region (out of the array bound), but somehow still belonging to the process. Thus, it does not throw segmentation fault and it will give garbage values instead.

To fix the issue, I recommend the following change (note -1 in rpad).

class MySource(Source):

    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.count = 1

    @property
    def count(self):
        """Return the number of sources/receives."""
        return self._count

    @count.setter
    def count(self, count):
        """Return the number of sources/receives."""
        self._count = count

    @property
    def interpolated_points_and_values(self):
        """ Calculate the amplitude value at all points"""

        points = np.array([], dtype=np.uint)
        values = np.array([], dtype=self.space_model.dtype)
        offset = np.array([0], dtype=np.uint)

        nbl_pad_width = self.space_model.nbl_pad_width
        halo_pad_width = self.space_model.halo_pad_width
        dimension = self.space_model.dimension
        bbox = self.space_model.bounding_box

        # senoid window
        for dim, dim_length in enumerate(self.space_model.shape):

            # points
            lpad = halo_pad_width[dim][0] + nbl_pad_width[dim][0]
            rpad = halo_pad_width[dim][1] + nbl_pad_width[dim][1] + dim_length - 1
            points = np.append(points, np.array([lpad, rpad], dtype=np.uint))
            # values
            coor_min, coor_max = bbox[dim * dimension:2 + dim * dimension]
            coordinates = np.linspace(coor_min, coor_max, dim_length)
            amplitudes = np.sin(np.pi / (coor_max - coor_min) * coordinates)
            values = np.append(values, amplitudes.astype(values.dtype))
        # offset
        offset = np.append(offset, np.array([values.size], dtype=np.uint))

        # add window for density term

        # points
        lpad_x = halo_pad_width[0][0] + nbl_pad_width[0][0]
        rpad_x = halo_pad_width[0][1] + nbl_pad_width[0][1] + self.space_model.shape[0] - 1
        points = np.append(points, np.array([lpad_x, rpad_x], dtype=np.uint))

        lpad_z = halo_pad_width[1][0] + nbl_pad_width[1][0]
        rpad_z = halo_pad_width[1][1] + nbl_pad_width[1][1] + self.space_model.shape[1] - 1
        points = np.append(points, np.array([lpad_z, rpad_z], dtype=np.uint))

        xmin, xmax = bbox[0: 2]
        zmin, zmax = bbox[2: 4]
        kx = np.pi / (xmax - xmin)
        kz = np.pi / (zmax - zmin)
        x_coordinates = np.linspace(xmin, xmax, self.space_model.shape[0])
        z_coordinates = np.linspace(zmin, zmax, self.space_model.shape[1])
        X, Z = np.meshgrid(z_coordinates, z_coordinates)

        # (1 / rho) * drho_dz * dphi_dz
        x_amplitudes = np.sin(kx * x_coordinates)
        drho_dz = np.gradient(dens[:, 256], z_coordinates).mean()
        z_amplitudes = kz * drho_dz * np.cos(kz * z_coordinates) / dens[:,256]

        values = np.append(values, x_amplitudes.astype(values.dtype))
        values = np.append(values, z_amplitudes.astype(values.dtype))

        offset = np.append(offset, np.array([values.size], dtype=np.uint))

        return points, values, offset

[joao-bapdm]

Thanks @jaimesouza, I've made the changes.