Use the convenient formulation for 3d and variable density case too (#60)

Author: joao-bapdm

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

@@ -177,8 +177,8 @@ double forward(f_type *u, f_type *velocity, f_type *damp,
                     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 numerator = (1.0 - damp[domain_offset] * dt / (2 * slowness));
+                    f_type denominator = (1.0 + damp[domain_offset] * dt / 2 );
+                    f_type numerator = (1.0 - damp[domain_offset] * dt / 2 );
 
                     value *= (dtSquared / slowness) / denominator;
 
@@ -272,7 +272,7 @@ double forward(f_type *u, f_type *velocity, f_type *damp,
                             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 denominator = (1.0 + damp[domain_offset] * dt / 2 );
 
                             f_type value = dtSquared / slowness * kws * wavelet[wavelet_offset] / denominator;
 

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

@@ -192,8 +192,8 @@ double forward(f_type *u, f_type *velocity, f_type *density, f_type *damp,
                 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 numerator = (1.0 - damp[domain_offset] * dt / (2 * slowness));
+                f_type denominator = (1.0 + damp[domain_offset] * dt / 2 );
+                f_type numerator = (1.0 - damp[domain_offset] * dt / 2 );
 
                 value *= (dtSquared / slowness) / denominator;
 
@@ -273,7 +273,7 @@ double forward(f_type *u, f_type *velocity, f_type *density, f_type *damp,
                         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 denominator = (1.0 + damp[domain_offset] * dt / 2 );
 
                         f_type value = dtSquared / slowness * kws * wavelet[wavelet_offset] / denominator;
 

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

@@ -203,8 +203,8 @@ double forward(f_type *u, f_type *velocity, f_type *density, f_type *damp,
                     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 numerator = (1.0 - damp[domain_offset] * dt / (2 * slowness));
+                    f_type denominator = (1.0 + damp[domain_offset] * dt / 2 );
+                    f_type numerator = (1.0 - damp[domain_offset] * dt / 2 );
 
                     value *= (dtSquared / slowness) / denominator;
 
@@ -297,7 +297,7 @@ double forward(f_type *u, f_type *velocity, f_type *density, f_type *damp,
                             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 denominator = (1.0 + damp[domain_offset] * dt / 2 );
 
                             f_type value = dtSquared / slowness * kws * wavelet[wavelet_offset] / denominator;
 

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

@@ -151,8 +151,6 @@ double gradient(f_type *u, f_type *v, f_type *grad, f_type *velocity, f_type *da
                 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 numerator = (1.0 - damp[domain_offset] * dt / (2 * slowness));
                 f_type denominator = (1.0 + damp[domain_offset] * dt / 2 );
                 f_type numerator = (1.0 - damp[domain_offset] * dt / 2 );
 
@@ -234,7 +232,6 @@ double gradient(f_type *u, f_type *v, f_type *grad, f_type *velocity, f_type *da
                         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 denominator = (1.0 + damp[domain_offset] * dt / 2 );
 
                         f_type value = dtSquared / slowness * kws * wavelet[wavelet_offset] / denominator;

tests/reference_solution/wavefield-3d-so-2.npy

@@ -32,7 +32,7 @@ def test_properties(self, dimension):
         assert space_model.damping_length == (0.0,) * dimension * 2
         assert space_model.boundary_condition == ('none',) * dimension * 2
         assert space_model.damping_polynomial_degree == 3
-        assert space_model.damping_alpha == 0.001
+        assert space_model.damping_alpha == 0.002
 
     @pytest.mark.parametrize(
         'dimension, bbox, spacing, shape', [