Examples: reworked 2nd order acoustic example to use new stencil evaluation and API

Author: EdCaunt

devitoboundary/stencil_cache.dat

@@ -282,8 +282,9 @@ def evaluate_stencils(df, point_type, n_stencils, left_variants, right_variants,
 
     if point_type == 'paired_right':
         dst = df.dist.to_numpy()[:, np.newaxis]
+        # dst used to have a minus (this was wrong)
         eta_left = np.tile(df.eta_l.to_numpy()[:, np.newaxis],
-                           (1, n_stencils)) + n_stencils - dst - eta_base - 1
+                           (1, n_stencils)) + n_stencils + dst - eta_base - 1
         eta_right = np.tile(df.eta_r.to_numpy()[:, np.newaxis],
                             (1, n_stencils)) + n_stencils - eta_base - 1
         for left_var in range(space_order+1):
@@ -294,7 +295,6 @@ def evaluate_stencils(df, point_type, n_stencils, left_variants, right_variants,
                     func = coeff_functions[left_var, right_var, coeff]
                     stencils[mask, coeff] = func(eta_left[mask],
                                                  eta_right[mask])
-
     return stencils
 
 
@@ -531,11 +531,10 @@ def get_component_weights(data, axis, function, deriv, stencil_generator):
     w_shape = f_grid.shape + (ncoeffs,)
     w_dims = f_grid.dimensions + (s_dim,)
 
-    w = Function(name='w', dimensions=w_dims, shape=w_shape)
+    w = Function(name='w_'+axis_dim, dimensions=w_dims, shape=w_shape)
 
     w.data[:] = standard_stencil(deriv, function.space_order)
 
-    # FIXME: All of these except "double" have a bug
     # Fill the stencils
     get_variants(first, function.space_order, 'first',
                  axis_dim, stencil_generator, w)
@@ -553,7 +552,7 @@ def get_component_weights(data, axis, function, deriv, stencil_generator):
         get_variants(paired_right, function.space_order, 'paired_right',
                      axis_dim, stencil_generator, w)
 
-    w.data[:] /= f_grid.spacing[axis]  # Divide everything through by spacing
+    w.data[:] /= f_grid.spacing[axis]**deriv  # Divide everything through by spacing
 
     return w
 
@@ -589,13 +588,15 @@ def get_weights(data, function, deriv, bcs, offsets=(0, 0, 0)):
                           stencil_file=stencil_file)
 
     # FIXME: This will want to cope with varying numbers of dims in the future
+    # FIXME: Why does the name fix work?
     weights = [None for i in range(3)]
     for axis in range(3):
         sten_gen.all_variants(deriv, offsets[axis])
         axis_weights = get_component_weights(data[axis].data, axis, function,
                                              deriv, sten_gen)
+        # Am I doing something dumb here?
+        print(deriv, function, function.grid.dimensions[axis], axis_weights)
         weights[axis] = Coefficient(deriv, function,
                                     function.grid.dimensions[axis],
                                     axis_weights)
-
     return Substitutions(*tuple(weights))

devitoboundary/stencils/evaluation.py

@@ -83,7 +83,7 @@ class ImmersedBoundary:
     ----------
     name : str
         The name of the boundary surface.
-    suface : str
+    surface : str
         The path to the geometry file
     functions : pandas DataFrame
         A dataframe of the functions to which the immersed boundary surface is

devitoboundary/topography.py

@@ -1,5 +1,5 @@
-import matplotlib.pyplot as plt
 import numpy as np
+import pandas as pd
 
 from devito import Grid, TimeFunction, Eq, solve, Operator, ConditionalDimension
 from devitoboundary import ImmersedBoundary
@@ -42,12 +42,14 @@
 # Surface configuration
 infile = 'topography/crater_lake.ply'
 # Zero even derivatives on the boundary
-bc_0 = Eq(generic_function(x_b), 0)  # u(x_b) = 0
-bc_2 = Eq(generic_function(x_b, 2), 0)  # d^2u/dx^2(x_b) = 0
-bcs = {u: [bc_0, bc_2]}
+bcs_u = [Eq(generic_function(x_b, 2*i), 0)
+         for i in range(1+u.space_order//2)]
+functions = pd.DataFrame({'function': [u],
+                          'bcs': [bcs_u]},
+                         columns=['function', 'bcs'])
 
 # Create the immersed boundary surface
-surface = ImmersedBoundary(infile, u, bcs)
+surface = ImmersedBoundary('topography', infile, functions)
 
 # Configure the source
 time_range = TimeAxis(start=t0, stop=tn, step=dt)
@@ -60,11 +62,14 @@
 src.coordinates.data[0, -1] = -500  # 500m below sea level
 
 # Configure derivative needed
-deriv = ('u.d2',)
+derivs = pd.DataFrame({'function': [u],
+                       'derivative': [2]},
+                      columns=['function', 'derivative'])
+coeffs = surface.subs(derivs)
 
 # We can now write the PDE
 pde = VP*u.dt2 - u.laplace
-eq = Eq(pde, 0, coefficients=surface.subs(deriv))
+eq = Eq(pde, 0, coefficients=coeffs['substitution'].values[0])
 
 # And set up the update
 stencil = solve(eq.evaluate, u.forward)
@@ -78,9 +83,9 @@
 # And run
 op.apply(dt=dt)
 
-# outfile = 'data/seismic_topography_wavefield.npy'
-# np.save(outfile, usave.data)
-
+outfile = 'data/seismic_topography_wavefield.npy'
+np.save(outfile, usave.data)
+"""
 plot_extent = [0, grid.extent[0],
                origin[2], grid.extent[2] + origin[2]]
 for i in range(usave.data.shape[0] - 1):
@@ -94,3 +99,4 @@
     plt.savefig("images/image-%s" % str(i))
     # plt.show()
     plt.close()
+"""

examples/seismic_topography_example.py

@@ -65,13 +65,20 @@ def test_type_splitting(self, axis):
         assert(np.all(paired_left.index.get_level_values(xyz[axis]).to_numpy() == 3))
         assert(np.all(paired_right.index.get_level_values(xyz[axis]).to_numpy() == 5))
 
+
+class TestStencils:
+    """
+    A class containing tests to check stencil evaluation.
+    """
+    # FIXME: Need to check evaluate_stencils
+    # FIXME: Need to check fill_weights
+    # FIXME: Wants to check several grid spacings
     @pytest.mark.parametrize('axis', [0, 1, 2])
     @pytest.mark.parametrize('deriv', [1, 2])
     def test_stencil_evaluation(self, axis, deriv):
         """
         A test to check that stencils are evaluated to their correct values.
         """
-        # FIXME: In wrong class
         def evaluate_variant(stencil_generator, left_var, right_var,
                              left_eta, right_eta):
             """Evaluate the specified stencil"""
@@ -132,7 +139,7 @@ def check_row(data, index, stencil_generator):
 
         sten_gen = StencilGen(function.space_order, bcs,
                               stencil_file=stencil_file)
-        sten_gen.all_variants(deriv, 0)
+        sten_gen.all_variants(deriv, 0.)
 
         w = get_component_weights(distances, axis, function, deriv, sten_gen)