How to use the devito.Eq function in devito

"""
        Check the shape of the Array used to store an aliasing expression.
        The shape is impacted by loop blocking, which reduces the required
        write-to space.
        """
        grid = Grid(shape=(3, 3, 3))
        x, y, z = grid.dimensions  # noqa
        t = grid.stepping_dim

        f = Function(name='f', grid=grid)
        f.data_with_halo[:] = 1.
        u = TimeFunction(name='u', grid=grid, space_order=3)
        u.data_with_halo[:] = 0.5

        # Leads to 3D aliases
        eqn = Eq(u.forward, ((u[t, x, y, z] + u[t, x+1, y+1, z+1])*3*f +
                             (u[t, x+2, y+2, z+2] + u[t, x+3, y+3, z+3])*3*f + 1))
        op0 = Operator(eqn, opt=('noop', {'openmp': True}))
        op1 = Operator(eqn, opt=('advanced', {'openmp': True, 'cire-mincost-sops': 1}))

        x0_blk_size = op1.parameters[2]
        y0_blk_size = op1.parameters[3]
        z_size = op1.parameters[4]

        # Check Array shape
        arrays = [i for i in FindSymbols().visit(op1._func_table['bf0'].root)
                  if i.is_Array and i._mem_local]
        assert len(arrays) == 1
        a = arrays[0]
        assert len(a.dimensions) == 3
        assert a.halo == ((1, 1), (1, 1), (1, 1))
        assert Add(*a.symbolic_shape[0].args) == x0_blk_size + 2

def test_staggered(self, ndim):
        """
        Test the "deformed" allocation for staggered functions
        """
        grid = Grid(shape=tuple([11]*ndim))
        for stagg in tuple(powerset(grid.dimensions))[1::] + (NODE, CELL):
            f = Function(name='f', grid=grid, staggered=stagg)
            assert f.data.shape == tuple([11]*ndim)
            # Add a non-staggered field to ensure that the auto-derived
            # dimension size arguments are at maximum
            g = Function(name='g', grid=grid)
            # Test insertion into a central point
            index = tuple(5 for _ in f.dimensions)
            set_f = Eq(f[index], 2.)
            set_g = Eq(g[index], 3.)

            Operator([set_f, set_g])()
            assert f.data[index] == 2.

loop nests respectively.
        """
        grid = Grid(shape=(3, 3, 3))
        x, y, z = grid.dimensions  # noqa
        t = grid.stepping_dim
        i = Dimension(name='i')

        f = Function(name='f', grid=grid)
        f.data_with_halo[:] = 1.
        g = Function(name='g', shape=(3,), dimensions=(i,))
        g.data[:] = 2.
        u = TimeFunction(name='u', grid=grid, space_order=3)
        v = TimeFunction(name='v', grid=grid, space_order=3)

        # Leads to 3D aliases
        eqns = [Eq(u.forward, ((u[t, x, y, z] + u[t, x+1, y+1, z+1])*3*f +
                               (u[t, x+2, y+2, z+2] + u[t, x+3, y+3, z+3])*3*f + 1)),
                Inc(u[t+1, i, i, i], g + 1),
                Eq(v.forward, ((v[t, x, y, z] + v[t, x+1, y+1, z+1])*3*u.forward +
                               (v[t, x+2, y+2, z+2] + v[t, x+3, y+3, z+3])*3*u.forward +
                               1))]
        op0 = Operator(eqns, opt=('noop', {'openmp': True}))
        op1 = Operator(eqns, opt=('advanced', {'openmp': True, 'cire-mincost-sops': 1}))

        # Check code generation
        assert 'bf0' in op1._func_table
        assert 'bf1' in op1._func_table
        trees = retrieve_iteration_tree(op1._func_table['bf0'].root)
        assert len(trees) == 2
        assert trees[0][-1].nodes[0].body[0].write.is_Array
        assert trees[1][-1].nodes[0].body[0].write is u
        trees = retrieve_iteration_tree(op1._func_table['bf1'].root)

dims0 = Grid(shape).dimensions
        for dims in permutations(dims0):
            grid = Grid(shape=shape, dimensions=dims, dtype=np.float64)
            time = grid.time_dim
            a = TimeFunction(name='a', grid=grid, time_order=2, space_order=2)
            p_aux = Dimension(name='p_aux')
            b = Function(name='b', shape=shape + (10,), dimensions=dims + (p_aux,),
                         space_order=2, dtype=np.float64)
            b.data_with_halo[:] = 1.0
            b2 = Function(name='b2', shape=(10,) + shape, dimensions=(p_aux,) + dims,
                          space_order=2, dtype=np.float64)
            b2.data_with_halo[:] = 1.0
            eqns = [Eq(a.forward, a.laplace + 1.),
                    Eq(b, time*b*a + b)]
            eqns2 = [Eq(a.forward, a.laplace + 1.),
                     Eq(b2, time*b2*a + b2)]
            subs = {d.spacing: v for d, v in zip(dims0, [2.5, 1.5, 2.0][:grid.dim])}

            op = Operator(eqns, subs=subs, dle='noop')
            trees = retrieve_iteration_tree(op)
            assert len(trees) == 2
            assert all(trees[0][i] is trees[1][i] for i in range(3))

            op2 = Operator(eqns2, subs=subs, dle='noop')
            trees = retrieve_iteration_tree(op2)
            assert len(trees) == 2

            # Verify both operators produce the same result
            op(time=10)
            a.data_with_halo[:] = 0.
            op2(time=10)

bounds_ym[j] = j
            bounds_yM[j] = 2*n_domains-1-j

        bounds = (bounds_xm, bounds_xM, bounds_ym, bounds_yM)

        inner_sd = Inner(N=n_domains, bounds=bounds)

        grid = Grid(extent=(10, 10), shape=(10, 10), subdomains=(inner_sd, ))

        assert(grid.subdomains['inner'].shape == ((10, 1), (8, 1), (6, 1),
                                                  (4, 1), (2, 1)))

        f = TimeFunction(name='f', grid=grid, dtype=np.int32)
        f.data[:] = 0

        stencil = Eq(f.forward, solve(Eq(f.dt, 1), f.forward),
                     subdomain=grid.subdomains['inner'])

        op = Operator(stencil)
        op(time_m=0, time_M=9, dt=1)
        result = f.data[0]

        fex = Function(name='fex', grid=grid)
        expected = np.zeros((10, 10), dtype=np.int32)
        for j in range(0, n_domains):
            expected[j, j:10-j] = 10
        fex.data[:] = np.transpose(expected)

        assert((np.array(result) == np.array(fex.data[:])).all())

def test_constant_time_dense(self):
        """Test arithmetic between different data objects, namely Constant
        and Function."""
        i, j = dimify('i j')
        const = Constant(name='truc', value=2.)
        a = Function(name='a', shape=(20, 20), dimensions=(i, j))
        a.data[:] = 2.
        eqn = Eq(a, a + 2.*const)
        op = Operator(eqn)

        op.apply(a=a, truc=const)
        assert(np.allclose(a.data, 6.))

        # Applying a different constant still works
        op.apply(a=a, truc=Constant(name='truc2', value=3.))
        assert(np.allclose(a.data, 12.))

def abc_eq(self, abc_dim, left=True):
        """
        Equation of the absorbing boundary condition as a complement of the PDE
        :param val: symbolic value of the dampening profile
        :return: Symbolic equation inside the boundary layer
        """
        # Define time sep to be updated
        next = self.field.forward if self.forward else self.field.backward
        # Define PDE
        eta = self.damp_profile_init(abc_dim, left=left)
        eq = self.m * self.field.dt2 - self.field.laplace
        eq += eta * self.field.dt if self.forward else -eta * self.field.dt
        # Solve the symbolic equation for the field to be updated
        eq_time = solve(eq, next, rational=False, simplify=False)[0]
        # return the Stencil with H replaced by its symbolic expression
        return Eq(next, eq_time).subs({abc_dim.parent: abc_dim})

# Initalise
u.data[:] = 0.0
v.data[:] = 0.0

# Modified coefficients
#u_x_coeffs = (1, u, x[0], np.array([1.0, -2.0, 1.0]))
#u_t_coeffs = (1, u, time, np.array([1.0, -2.0, 1.0]))

u_x_coeffs = (1, u, x[0], np.array([-0.5, 0.0, 0.5]))
u_t_coeffs = (1, u, time, np.array([-0.5, 0.0, 0.5]))

coeffs=Coefficients(u_x_coeffs,u_t_coeffs)
#coeffs=Coefficients(u_x_coeffs)

# Main equations
eq = Eq(u.dt+u.dx+v.dx, coefficients=coeffs)
#eq = Eq(u.dt+u.dx+v.dx)

print(eq)
#help(eq)

def g3_tilde(field):
    return field.dz(x0=z-z.spacing/2)


# Functions  for additional factorization
b1mf = Function(name='b1mf', grid=grid, space_order=space_order)
b1m2e = Function(name='b1m2e', grid=grid, space_order=space_order)
b1mfa2 = Function(name='b1mfa2', grid=grid, space_order=space_order)
b1mfpfa2 = Function(name='b1mfpfa2', grid=grid, space_order=space_order)
bfes1ma2 = Function(name='bfes1ma2', grid=grid, space_order=space_order)

# Equations for additional factorization
eq_b1mf = Eq(b1mf, b * (1 - f))
eq_b1m2e = Eq(b1m2e, b * (1 + 2 * eps))
eq_b1mfa2 = Eq(b1mfa2, b * (1 - f * eta**2))
eq_b1mfpfa2 = Eq(b1mfpfa2, b * (1 - f + f * eta**2))
eq_bfes1ma2 = Eq(bfes1ma2, b * f * eta * sqrt(1 - eta**2))

# Time update equation for quasi-P state variable p
update_p_nl = t.spacing**2 * vel**2 / b * \
    (g1_tilde(b1m2e * g1(p0)) +
     g2_tilde(b1m2e * g2(p0)) +
     g3_tilde(b1mfa2 * g3(p0) + bfes1ma2 * g3(m0))) + \
    (2 - t.spacing * wOverQ) * p0 + (t.spacing * wOverQ - 1) * p0.backward

# Time update equation for quasi-S state variable m
update_m_nl = t.spacing**2 * vel**2 / b * \
    (g1_tilde(b1mf * g1(m0)) +
     g2_tilde(b1mf * g2(m0)) +
     g3_tilde(b1mfpfa2 * g3(m0) + bfes1ma2 * g3(p0))) + \
    (2 - t.spacing * wOverQ) * m0 + (t.spacing * wOverQ - 1) * m0.backward

b1mfa2 = Function(name='b1mfa2', grid=grid, space_order=space_order)
b1mfpfa2 = Function(name='b1mfpfa2', grid=grid, space_order=space_order)
bfes1ma2 = Function(name='bfes1ma2', grid=grid, space_order=space_order)

p_x = Function(name='p_x', grid=grid, space_order=space_order)
p_y = Function(name='p_y', grid=grid, space_order=space_order)
p_z = Function(name='p_z', grid=grid, space_order=space_order)
m_x = Function(name='m_x', grid=grid, space_order=space_order)
m_y = Function(name='m_y', grid=grid, space_order=space_order)
m_z = Function(name='m_z', grid=grid, space_order=space_order)

# Equations for additional factorization
eq_b1mf = Eq(b1mf, b * (1 - f))
eq_b1m2e = Eq(b1m2e, b * (1 + 2 * eps))
eq_b1mfa2 = Eq(b1mfa2, b * (1 - f * eta**2))
eq_b1mfpfa2 = Eq(b1mfpfa2, b * (1 - f + f * eta**2))
eq_bfes1ma2 = Eq(bfes1ma2, b * f * eta * sqrt(1 - eta**2))

eq_px = Eq(p_x, g1_tilde(b1m2e * g1(p0)))
eq_py = Eq(p_y, g2_tilde(b1m2e * g2(p0)))
eq_pz = Eq(p_z, g3_tilde(b1mfa2 * g3(p0) + bfes1ma2 * g3(m0)))

eq_mx = Eq(m_x, g1_tilde(b1mf * g1(m0)))
eq_my = Eq(m_y, g2_tilde(b1mf * g2(m0)))
eq_mz = Eq(m_z, g3_tilde(b1mfpfa2 * g3(m0) + bfes1ma2 * g3(p0)))

# Time update equation for quasi-P state variable p
update_p_nl = t.spacing**2 * vel**2 / b * (p_x + p_y + p_z) + \
    (2 - t.spacing * wOverQ) * p0 + (t.spacing * wOverQ - 1) * p0.backward

# Time update equation for quasi-S state variable m
update_m_nl = t.spacing**2 * vel**2 / b * (m_x + m_y + m_z) + \