Set initial conditions

[joao-bapdm]

So that the initial pressure, or initial "pressure velocity" can be non-zero

[jaimesouza]

I don’t understand. Where should it happen?

[krober10nd]

probably in u_full, recv = solver.forward()

you could pass a u0 keyword which sets the initial pressure field.

[krober10nd]

like solver.forward(u0=initial_u)

[jaimesouza]

I see. But why would we need that right now?

[krober10nd]

in order to validate the numerical implementation with the method of manufactured solutions one needs to set an initial condition for the pressure field.

[krober10nd]

analytical solutions exist for the constant density but not for the variable one to my knowledge. the MMS allows us to validate all schemes.

[joao-bapdm]

I see. But why would we need that right now?

We don't need it right now. I would have put an "enhancement" label if I could, or knew how to do it. The validation mentioned by @krober10nd is why this came to my attention however.

For now I'll try my luck with a manufactured solution with homogeneous initial conditions.

[krober10nd]

oh that won't work. The initial condition is the forcing mechanism.

[joao-bapdm]

I undestand there are two different ways to do a MMS:

• Set initial conditions but enforce no source, or forcing function of any kind
• Use a specific ansatz that results in a source that preserves homogeneous initial conditions. For instance, if the solution is some f(x) * g(t), and both g(t) and g'(t) are zero at t=0, then I can still use a field initially at rest.

[joao-bapdm]

Since u(x, 0) = f(x) g(0) = 0, and du/dt(x, 0) = f(x) dg/dt(0) = 0, I still have homogeneous initial conditions.

[krober10nd]

Okay, I've only done the first way. Lets see what happens.

[jaimesouza]

So, u0 is just a 2D/3D grid with the same shape as the extended u (only one timestep)? By extendend shape I mean (original + nbl).

[jaimesouza]

Or I just need original?

[joao-bapdm]

I think it should be in the extended grid.

[jaimesouza]

Right. It is on the queue of enhancements.

[joao-bapdm]

I am trying the initial conditions at add-initial-conditions. Here are some slight modifications to the 2D acoustic example.

from simwave import (
    SpaceModel, TimeModel, RickerWavelet, Solver, Compiler,
    Receiver, Source, plot_wavefield, plot_shotrecord, plot_velocity_model
)
import numpy as np
import matplotlib.pyplot as plt

# 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'
    # cflags='-O3 -fPIC -ffast-math -fopenmp \
    #       -fopenmp-targets=nvptx64-nvidia-cuda -Xopenmp-target -march=sm_75'
)

# 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=1
)

# create the set of sources
source = Source(
    space_model,
    coordinates=[(2560, 2560)],
    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, amplitude=0)

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

print("Timesteps:", time_model.timesteps)

# initial grid
initial_grid = np.zeros(shape=space_model.shape)
initial_grid[256,256] = 100

# run the forward
u_full, recv = solver.forward(initial_grid)

print("u_full shape:", u_full.shape)
plot_velocity_model(space_model.velocity_model)
plot_wavefield(u_full[-1])
plot_shotrecord(recv)

plt.plot(time_model.time_values, recv[:,256])
plt.title("Signal at $256^{th}$ receiver")
plt.xlabel("time [s]")
plt.ylabel("amplitude", rotation=90)
plt.show()

I've just set the ricker amplitude to zero, and added an initial value of 100 at the node (256,256) using initial_grid. I am asking for u_full back with saving_stride of 1, and looking at the initial values of it, I get:

u_full[0,:,:].max()
0.0

u_full[1,:,:].max()
0.0

u_full[2,:,:].max()
3.125

which doesn't seem right to me, since then initial_grid is not an actual initial condition.

[jaimesouza]

fixed in https://github.com/HPCSys-Lab/simwave/tree/add-initial-condition