Fix damping and examples (#70)

* damping

* damping

* computing eta

* alternative calculationa

* fixing forward

* fix model

* clean up model

* clean up model

* fix examples

* fix acoustic2d

Author: hermes-senger

examples/acoustic_2D.py

@@ -45,10 +45,11 @@
 )
 
 # Velocity model
+#vel = np.zeros(shape=(512, 512), dtype=np.float32)
 vel = np.zeros(shape=(512, 512), dtype=np.float32)
 vel[:] = 1500.0
-vel[100:] = 2000.0
-
+#vel[100:] = 2000.0
+#vel[1,1]=4500
 # create the space model
 space_model = SpaceModel(
     bounding_box=(0, 5120, 0, 5120),
@@ -61,14 +62,13 @@
 # config boundary conditions
 # (none,  null_dirichlet or null_neumann)
 space_model.config_boundary(
-    damping_length=(0, 510, 510, 510),
+        damping_length=(510, 510, 510, 510),
     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(

examples/acoustic_2D_density.py

@@ -66,13 +66,11 @@
 # config boundary conditions
 # (none,  null_dirichlet or null_neumann)
 space_model.config_boundary(
-    damping_length=0,
+    damping_length=(510, 510, 510, 510),
     boundary_condition=(
         "null_neumann", "null_dirichlet",
         "none", "null_dirichlet"
-    ),
-    damping_polynomial_degree=3,
-    damping_alpha=0.001
+    )
 )
 
 # create the time model

examples/acoustic_3D.py

@@ -59,16 +59,13 @@
 # config boundary conditions
 # (none,  null_dirichlet or null_neumann)
 space_model.config_boundary(
-    damping_length=0,
+    damping_length=(0, 100, 100, 100, 100, 100),
     boundary_condition=(
         "null_neumann", "null_dirichlet",
         "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,

examples/acoustic_3D_density.py

@@ -65,14 +65,12 @@
 # config boundary conditions
 # (none,  null_dirichlet or null_neumann)
 space_model.config_boundary(
-    damping_length=0,
+    damping_length=(0, 100, 100, 100, 100, 100),
     boundary_condition=(
         "null_neumann", "null_dirichlet",
         "null_dirichlet", "null_dirichlet",
         "null_dirichlet", "null_dirichlet"
-    ),
-    damping_polynomial_degree=3,
-    damping_alpha=0.001
+    )
 )
 
 # create the time model

examples/micro.py

@@ -0,0 +1,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)

setup.py

@@ -1,4 +1,3 @@
-import versioneer
 from setuptools import setup, find_packages
 
 with open("requirements.txt") as f:
@@ -9,8 +8,6 @@
 
 setup(
     name="simwave",
-    version=versioneer.get_version(),
-    cmdclass=versioneer.get_cmdclass(),
     description="Finite difference 2D/3D acoustic wave propagator.",
     long_description=readme,
     long_description_content_type="text/markdown",

simwave/kernel/backend/c_code/forward/constant_density/2d/wave.c

@@ -34,9 +34,11 @@ double forward(f_type *u, f_type *velocity, f_type *damp,
                size_t saving_stride, f_type dt,
                size_t begin_timestep, size_t end_timestep,
                size_t space_order, size_t num_snapshots){
-
+ 
     size_t stencil_radius = space_order / 2;
 
+    //size_t stencil_radius = space_order / 2;
+
     size_t domain_size = nz * nx;
 
     f_type dzSquared = dz * dz;
@@ -167,13 +169,16 @@ double forward(f_type *u, f_type *velocity, f_type *damp,
                 value += sum_x/dxSquared + sum_z/dzSquared;
 
                 // parameter to be used
-                f_type slowness = 1.0 / (velocity[domain_offset] * velocity[domain_offset]);
+                //f_type slowness = 1.0 / (velocity[domain_offset] * velocity[domain_offset]);
+                f_type c2 = velocity[domain_offset] * velocity[domain_offset];
 
                 // denominator with damp coefficient
-                f_type denominator = (1.0 + damp[domain_offset] * dt / (2 * slowness));
-                f_type numerator = (1.0 - damp[domain_offset] * dt / (2 * slowness));
+                // f_type denominator = (1.0 + damp[domain_offset] * dt / (2 * slowness));
+                f_type denominator = 1.0 + damp[domain_offset] * dt;
+                f_type numerator = 1.0 - damp[domain_offset] * dt;
 
-                value *= (dtSquared / slowness) / denominator;
+                //value *= (dtSquared / slowness) / denominator;
+                value *= (dtSquared * c2) / denominator;
 
                 u[next_snapshot] = 2.0 / denominator * u[current_snapshot] - (numerator / denominator) * u[prev_snapshot] + value;
             }

simwave/kernel/backend/c_code/forward/constant_density/3d/wave.c

@@ -174,13 +174,17 @@ double forward(f_type *u, f_type *velocity, f_type *damp,
                     value += sum_y/dySquared + sum_x/dxSquared + sum_z/dzSquared;
 
                     // parameter to be used
-                    f_type slowness = 1.0 / (velocity[domain_offset] * velocity[domain_offset]);
+                    //f_type slowness = 1.0 / (velocity[domain_offset] * velocity[domain_offset]);
+                    f_type c2 = velocity[domain_offset] * velocity[domain_offset];
 
                     // denominator with damp coefficient
-                    f_type denominator = (1.0 + damp[domain_offset] * dt / (2 * slowness));
-                    f_type numerator = (1.0 - damp[domain_offset] * dt / (2 * slowness));
+                    //f_type denominator = (1.0 + damp[domain_offset] * dt / (2 * slowness));
+                    f_type denominator = 1.0 + damp[domain_offset] * dt;
+                    //f_type numerator = (1.0 - damp[domain_offset] * dt / (2 * slowness));
+                    f_type numerator = 1.0 - damp[domain_offset] * dt;
 
-                    value *= (dtSquared / slowness) / denominator;
+                    //value *= (dtSquared /slowness) / denominator;
+                    value *= (dtSquared * c2) / denominator;
 
                     u[next_snapshot] = 2.0 / denominator * u[current_snapshot] - (numerator / denominator) * u[prev_snapshot] + value;
                 }

simwave/kernel/backend/c_code/forward/variable_density/2d/wave.c

@@ -189,13 +189,17 @@ double forward(f_type *u, f_type *velocity, f_type *density, f_type *damp,
                 value -= (term_x + term_z) / density[domain_offset];
 
                 // parameter to be used
-                f_type slowness = 1.0 / (velocity[domain_offset] * velocity[domain_offset]);
+                //f_type slowness = 1.0 / (velocity[domain_offset] * velocity[domain_offset]);
+                f_type c2 = velocity[domain_offset] * velocity[domain_offset];
 
                 // denominator with damp coefficient
-                f_type denominator = (1.0 + damp[domain_offset] * dt / (2 * slowness));
-                f_type numerator = (1.0 - damp[domain_offset] * dt / (2 * slowness));
+                //f_type denominator = (1.0 + damp[domain_offset] * dt / (2 * slowness));
+                f_type denominator = 1.0 + damp[domain_offset] * dt;
+                //f_type numerator = (1.0 - damp[domain_offset] * dt / (2 * slowness));
+                f_type numerator = 1.0 - damp[domain_offset] * dt;
 
-                value *= (dtSquared / slowness) / denominator;
+                //value *= (dtSquared / slowness) / denominator;
+                value *= (dtSquared * c2) / denominator;
 
                 u[next_snapshot] = 2.0 / denominator * u[current_snapshot] - (numerator / denominator) * u[prev_snapshot] + value;
             }

simwave/kernel/backend/c_code/forward/variable_density/3d/wave.c

@@ -200,13 +200,17 @@ double forward(f_type *u, f_type *velocity, f_type *density, f_type *damp,
                     value -= (term_y + term_x + term_z) / density[domain_offset];
 
                     // parameter to be used
-                    f_type slowness = 1.0 / (velocity[domain_offset] * velocity[domain_offset]);
+                    //f_type slowness = 1.0 / (velocity[domain_offset] * velocity[domain_offset]);
+                    f_type c2 = velocity[domain_offset] * velocity[domain_offset];
 
                     // denominator with damp coefficient
-                    f_type denominator = (1.0 + damp[domain_offset] * dt / (2 * slowness));
-                    f_type numerator = (1.0 - damp[domain_offset] * dt / (2 * slowness));
+		    // f_type denominator = (1.0 + damp[domain_offset] * dt / (2 * slowness));
+                    f_type denominator = 1.0 + damp[domain_offset] * dt;
+                    // f_type numerator = (1.0 - damp[domain_offset] * dt / (2 * slowness));
+                    f_type numerator = 1.0 - damp[domain_offset] * dt;
 
-                    value *= (dtSquared / slowness) / denominator;
+                    //value *= (dtSquared / slowness) / denominator;
+                    value *= (dtSquared * c2) / denominator;
 
                     u[next_snapshot] = 2.0 / denominator * u[current_snapshot] - (numerator / denominator) * u[prev_snapshot] + value;
                 }