simwave/examples/micro.py at d2dad4b02a04ad3e3ca2bd0a271d33d70573886f · HPCSys-Lab/simwave · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
from simwave import (
    SpaceModel, TimeModel, RickerWavelet, Solver, Compiler,
    Receiver, Source, plot_wavefield, plot_shotrecord, plot_velocity_model
)
import numpy as np


# available language options:
# c (sequential)
# cpu_openmp (parallel CPU)
# gpu_openmp (GPU)
# gpu_openacc (GPU)
compiler_options = {
    'c': {
        'cc': 'gcc',
        'language': 'c',
        'cflags': '-O3 -fPIC -ffast-math -Wall -std=c99 -shared'
    },
    'cpu_openmp': {
        'cc': 'gcc',
        'language': 'cpu_openmp',
        'cflags': '-O3 -fPIC -ffast-math -Wall -std=c99 -shared -fopenmp'
    },
    'gpu_openmp': {
        'cc': 'clang',
        'language': 'gpu_openmp',
        'cflags': '-O3 -fPIC -ffast-math -fopenmp \
                   -fopenmp-targets=nvptx64-nvidia-cuda \
                   -Xopenmp-target -march=sm_75'
    },
    'gpu_openacc': {
        'cc': 'pgcc',
        'language': 'gpu_openacc',
        'cflags': '-O3 -fPIC -acc:gpu -gpu=pinned -mp'
    },
}

selected_compiler = compiler_options['c']

# set compiler options
compiler = Compiler(
    cc=selected_compiler['cc'],
    language=selected_compiler['language'],
    cflags=selected_compiler['cflags']
)

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

# create the space model
space_model = SpaceModel(
    bounding_box=(0, 20, 0, 20),
    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, 1010, 1010, 1010),
        damping_length=(20, 20, 20, 20),
#    damping_length=0,
    boundary_condition=(
        #        "null_neumann", "null_dirichlet",
        "null_dirichlet", "null_dirichlet",
        "null_dirichlet", "null_dirichlet"
    )
)

print(' damping_alpha=',space_model.damping_alpha)

# create the time model
time_model = TimeModel(
    space_model=space_model,
    tf=0.01,
    saving_stride=0
)

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

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

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

print("Timesteps:", time_model.timesteps)

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