Fix tests and implement the new formulation in variable density code

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

@@ -169,7 +169,7 @@ double forward(f_type *u, f_type *velocity, f_type *damp,
                 // parameter to be used
                 f_type slowness = 1.0 / (velocity[domain_offset] * velocity[domain_offset]);
 
-                //denominator with damp coefficient
+                // 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));
 
@@ -249,7 +249,7 @@ double forward(f_type *u, f_type *velocity, f_type *damp,
                         // parameter to be used
                         f_type slowness = 1.0 / (velocity[domain_offset] * velocity[domain_offset]);
 
-                        //denominator with damp coefficient
+                        // denominator with damp coefficient
                         f_type denominator = (1.0 + damp[domain_offset] * dt / (2 * slowness));
 
                         f_type value = dtSquared / slowness * kws * wavelet[wavelet_offset] / denominator;

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

@@ -176,7 +176,7 @@ double forward(f_type *u, f_type *velocity, f_type *damp,
                     // parameter to be used
                     f_type slowness = 1.0 / (velocity[domain_offset] * velocity[domain_offset]);
 
-                    //denominator with damp coefficient
+                    // 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));
 
@@ -271,7 +271,7 @@ double forward(f_type *u, f_type *velocity, f_type *damp,
                             // parameter to be used
                             f_type slowness = 1.0 / (velocity[domain_offset] * velocity[domain_offset]);
 
-                            //denominator with damp coefficient
+                            // denominator with damp coefficient
                             f_type denominator = (1.0 + damp[domain_offset] * dt / (2 * slowness));
 
                             f_type value = dtSquared / slowness * kws * wavelet[wavelet_offset] / denominator;

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

@@ -188,12 +188,16 @@ double forward(f_type *u, f_type *velocity, f_type *density, f_type *damp,
 
                 value -= (term_x + term_z) / density[domain_offset];
 
-                //denominator with damp coefficient
-                f_type denominator = (1.0 + damp[domain_offset] * dt);
+                // parameter to be used
+                f_type slowness = 1.0 / (velocity[domain_offset] * velocity[domain_offset]);
 
-                value *= (dtSquared * velocity[domain_offset] * velocity[domain_offset]) / denominator;
+                // 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));
 
-                u[next_snapshot] = 2.0 / denominator * u[current_snapshot] - ((1.0 - damp[domain_offset] * dt) / denominator) * u[prev_snapshot] + value;
+                value *= (dtSquared / slowness) / denominator;
+
+                u[next_snapshot] = 2.0 / denominator * u[current_snapshot] - (numerator / denominator) * u[prev_snapshot] + value;
             }
         }
 
@@ -265,7 +269,13 @@ double forward(f_type *u, f_type *velocity, f_type *density, f_type *damp,
                         size_t domain_offset = i * nx + j;
                         size_t next_snapshot = next_t * domain_size + domain_offset;
 
-                        f_type value = dtSquared * velocity[domain_offset] * velocity[domain_offset] * kws * wavelet[wavelet_offset];
+                        // parameter to be used
+                        f_type slowness = 1.0 / (velocity[domain_offset] * velocity[domain_offset]);
+
+                        // denominator with damp coefficient
+                        f_type denominator = (1.0 + damp[domain_offset] * dt / (2 * slowness));
+
+                        f_type value = dtSquared / slowness * kws * wavelet[wavelet_offset] / denominator;
 
                         #if defined(CPU_OPENMP) || defined(GPU_OPENMP)
                         #pragma omp atomic

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

@@ -199,12 +199,16 @@ double forward(f_type *u, f_type *velocity, f_type *density, f_type *damp,
 
                     value -= (term_y + term_x + term_z) / density[domain_offset];
 
-                    //denominator with damp coefficient
-                    f_type denominator = (1.0 + damp[domain_offset] * dt);
+                    // parameter to be used
+                    f_type slowness = 1.0 / (velocity[domain_offset] * velocity[domain_offset]);
 
-                    value *= (dtSquared * velocity[domain_offset] * velocity[domain_offset]) / denominator;
+                    // 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));
 
-                    u[next_snapshot] = 2.0 / denominator * u[current_snapshot] - ((1.0 - damp[domain_offset] * dt) / denominator) * u[prev_snapshot] + value;
+                    value *= (dtSquared / slowness) / denominator;
+
+                    u[next_snapshot] = 2.0 / denominator * u[current_snapshot] - (numerator / denominator) * u[prev_snapshot] + value;
                 }
             }
         }
@@ -287,9 +291,15 @@ double forward(f_type *u, f_type *velocity, f_type *density, f_type *damp,
 
                             // current source point in the grid
                             size_t domain_offset = (i * nx + j) * ny + k;
-                            size_t next_snapshot = next_t * domain_size + domain_offset;
+                            size_t next_snapshot = next_t * domain_size + domain_offset;                           
+
+                            // parameter to be used
+                            f_type slowness = 1.0 / (velocity[domain_offset] * velocity[domain_offset]);
+
+                            // denominator with damp coefficient
+                            f_type denominator = (1.0 + damp[domain_offset] * dt / (2 * slowness));
 
-                            f_type value = dtSquared * velocity[domain_offset] * velocity[domain_offset] * kws * wavelet[wavelet_offset];
+                            f_type value = dtSquared / slowness * kws * wavelet[wavelet_offset] / denominator;
 
                             #if defined(CPU_OPENMP) || defined(GPU_OPENMP)
                             #pragma omp atomic

simwave/kernel/frontend/source.py

@@ -72,8 +72,8 @@ def grid_positions(self):
                 zmin, zmax, xmin, xmax = self.space_model.bounding_box
                 z_spacing, x_spacing = self.space_model.grid_spacing
 
-                if not(zmin <= coord[0] <= zmax) or \
-                   not(xmin <= coord[1] <= xmax):
+                if not (zmin <= coord[0] <= zmax) or \
+                   not (xmin <= coord[1] <= xmax):
                     raise Exception("Coordinates %s out of bounds." % coord)
 
                 zpos = (coord[0] - zmin) / z_spacing
@@ -87,9 +87,9 @@ def grid_positions(self):
 
                 z_spacing, x_spacing, y_spacing = self.space_model.grid_spacing
 
-                if not(zmin <= coord[0] <= zmax) or \
-                   not(xmin <= coord[1] <= xmax) or \
-                   not(ymin <= coord[2] <= ymax):
+                if not (zmin <= coord[0] <= zmax) or \
+                   not (xmin <= coord[1] <= xmax) or \
+                   not (ymin <= coord[2] <= ymax):
                     raise Exception("Coordinates %s out of bounds." % coord)
 
                 zpos = (coord[0] - zmin) / z_spacing

tests/reference_solution/generator.py

@@ -1,5 +1,5 @@
-from pywave import SpaceModel, TimeModel, RickerWavelet, Solver
-from pywave import Receiver, Source
+from simwave import SpaceModel, TimeModel, RickerWavelet, Solver
+from simwave import Receiver, Source, Compiler
 import numpy as np
 
 
@@ -58,32 +58,38 @@ def generate_one(dimension, space_order):
     # create the time model
     time_model = TimeModel(
         space_model=space_model,
-        tf=tf
+        tf=tf,
+        saving_stride=0
     )
 
     # create the set of sources
     source = Source(
         space_model=space_model,
-        coordinates=source_position
+        coordinates=source_position,
+        window_radius=1
     )
 
     # crete the set of receivers
     receiver = Receiver(
         space_model=space_model,
-        coordinates=source_position
+        coordinates=source_position,
+        window_radius=1
     )
 
     # create a ricker wavelet with 10hz of peak frequency
     ricker = RickerWavelet(f0, time_model)
 
+    # create a compiler
+    compiler = Compiler(language='c')
+
     # create the solver
     solver = Solver(
+        compiler=compiler,
         space_model=space_model,
         time_model=time_model,
         sources=source,
         receivers=receiver,
-        wavelet=ricker,
-        saving_stride=0
+        wavelet=ricker
     )
 
     # run the forward

tests/test_solution.py

@@ -67,7 +67,7 @@ def test_solution(self, dimension, space_order, language, density):
 
         if density:
             den = np.zeros(shape=shape, dtype=np.float32)
-            den[:] = 5.0
+            den[:] = 1.0
 
             # create the space model
             space_model = SpaceModel(
@@ -136,4 +136,4 @@ def test_solution(self, dimension, space_order, language, density):
         # load the reference result
         u_ref = np.load(ref_file)
 
-        assert np.allclose(u, u_ref, atol=1e-04)
+        assert np.allclose(u, u_ref, atol=1e-05)