Symbolic FD coefficients fail with biharmonic

[ccuetom]

When we define a function to be symbolic,

u = TimeFunction(name="u", grid=model.grid, time_order=2, space_order=10, coefficients='symbolic')

it can be used on the second-order isotropic acoustic wave equation without problems,

pde = model.m * u.dt2 - u.laplace + model.damp * u.dt
Eq(u.forward, solve(pde, u.forward), coefficients=coeffs)

However, if the biharmonic is included, the operator cannot be generated,

pde = model.m * u.dt2 - u.laplace - dt**2/12 * u.biharmonic(1/model.m) + model.damp * u.dt
Eq(u.forward, solve(pde, u.forward), coefficients=coeffs)

[ccuetom]

Full MFE below:

import numpy as np
import sympy as sp
from devito import Grid, TimeFunction, Coefficient, Eq, Substitutions, solve,  Operator
from examples.seismic import RickerSource, TimeAxis
from examples.seismic import Model, plot_velocity


# Define symbol for laplacian replacement
H = sp.symbols('H')

# Define a physical size
Lx = 2000
Lz = Lx
h = 10
Nx = int(Lx/h)+1
Nz = Nx

shape = (Nx, Nz)  # Number of grid point
spacing = (h, h)  # Grid spacing in m. The domain size is now 2km by 2km
origin = (0., 0.)

# Define a velocity profile. The velocity is in km/s
v = np.empty(shape, dtype=np.float32)
v[:, :121] = 1.5
v[:, 121:] = 4.0

# With the velocity and model size defined, we can create the seismic model that
# encapsulates these properties. We also define the size of the absorbing layer as 10 grid points
nbl = 10
model = Model(vp=v, origin=origin, shape=shape, spacing=spacing,
              space_order=20, nbl=nbl, bcs="damp")

t0 = 0.  # Simulation starts a t=0
tn = 500.  # Simulation lasts 0.5 seconds (500 ms)
dt = 1.0  # Time step of 0.2ms

time_range = TimeAxis(start=t0, stop=tn, step=dt)

f0 = 0.015  # Source peak frequency is 25Hz (0.025 kHz)
src = RickerSource(name='src', grid=model.grid, f0=f0,
                   npoint=1, time_range=time_range)

# First, position source centrally in all dimensions, then set depth
src.coordinates.data[0, :] = np.array(model.domain_size) * .5
src.coordinates.data[0, -1] = 800.  # Depth is 800m

order = 10
u_DRP = TimeFunction(name="u_DRP", grid=model.grid, time_order=2, space_order=order, coefficients='symbolic')

pde_DRP = model.m * u_DRP.dt2 - H + model.damp * u_DRP.dt

# Define our custom FD coefficients:
x, z = model.grid.dimensions

# Lower layer
weights_l = np.array([  0.      ,  0.      ,  0.0274017,
                       -0.223818,  1.64875 , -2.90467,
                        1.64875 , -0.223818,  0.0274017,
                        0.      ,  0.       ])

ux_l_coeffs = Coefficient(2, u_DRP, x, weights_l/x.spacing**2)
uz_l_coeffs = Coefficient(2, u_DRP, z, weights_l/z.spacing**2)


# The underlying pde is the same in both subdomains
pde_DRP = model.m * u_DRP.dt2 - H + model.damp * u_DRP.dt

# Define our custom FD coefficients:
x, z = model.grid.dimensions
# Upper layer
weights_u = np.array([ 2.00462e-03, -1.63274e-02,  7.72781e-02,
                      -3.15476e-01,  1.77768e+00, -3.05033e+00,
                       1.77768e+00, -3.15476e-01,  7.72781e-02,
                      -1.63274e-02,  2.00462e-03])
# Lower layer
weights_l = np.array([  0.      ,  0.      ,  0.0274017,
                       -0.223818,  1.64875 , -2.90467,
                        1.64875 , -0.223818,  0.0274017,
                        0.      ,  0.       ])
# Create the Devito Coefficient objects:
ux_l_coeffs = Coefficient(2, u_DRP, x, weights_l/x.spacing**2)
uz_l_coeffs = Coefficient(2, u_DRP, z, weights_l/z.spacing**2)

# And the replacement rules:
coeffs_l = Substitutions(ux_l_coeffs,uz_l_coeffs)

# line below would work
laplace = u_DRP.laplace
# but this one wouldn't
laplace = u_DRP.laplace + dt**2/12 * u_DRP.biharmonic(1/model.m)

stencil_l = Eq(u_DRP.forward, solve(pde_DRP, u_DRP.forward).subs({H: laplace}), coefficients=coeffs_l)

# Source term:
src_term = src.inject(field=u_DRP.forward, expr=src * dt**2 / model.m)

# Create the operator
op = Operator([stencil_l] + src_term, subs=model.spacing_map)

op(time=time_range.num-1, dt=dt)

[FabioLuporini]

@ccuetom is this still an issue?

several symbolic coefficients tweaks have gone in recently

[ccuetom]

@FabioLuporini this seems to be partially solved: the biharmonic can now be computed, but the compilation fails with:

    op(time=time_range.num-1, dt=dt)
  File "devito/operator/operator.py", line 768, in __call__
    return self.apply(**kwargs)
  File "devito/operator/operator.py", line 839, in apply
    cfunction = self.cfunction
  File "devito/operator/operator.py", line 722, in cfunction
    self._lib = self._compiler.load(self._soname)
  File "devito/arch/compiler.py", line 260, in load
    return npct.load_library(str(self.get_jit_dir().joinpath(soname)), '.')
  File "lib/python3.10/site-packages/numpy/ctypeslib.py", line 158, in load_library
    return ctypes.cdll[libpath]
  File "lib/python3.10/ctypes/__init__.py", line 449, in __getitem__
    return getattr(self, name)
  File "lib/python3.10/ctypes/__init__.py", line 444, in __getattr__
    dll = self._dlltype(name)
  File "lib/python3.10/ctypes/__init__.py", line 374, in __init__
    self._handle = _dlopen(self._name, mode)
OSError: /tmp/devito-jitcache-uid1000/8a328b01328b5be5bdd1a0c4be261fa63e976661.so: undefined symbol: W

[mloubout]

@ccuetom you are using space_order=20 for the Model while using space_order=10 for u_DRP so it fails to find weights for the second part of the biharmonic because it is looking for 20-th order FD coefficients, if you switch the Model to space_order=10 to match the wavefield it works fine

[ccuetom]

@mloubout oh yeah, my bad!

[mloubout]

@ccuetom this is however a small bug as it should use the minimum order not maximum, fixed in

https://github.com/devitocodes/devito/pull/2254

that will allow you to use 20/10 orders (it will use 10 for derivatives using both but correctly)