How to use the devito.norm function in devito

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.

        eqn = Eq(u.forward, ((u[t, x, y] + u[t, x+2, y])*3*f +
                             (u[t, x+1, y+1] + u[t, x+3, y+1])*3*f + 1))
        op0 = Operator(eqn, dse='noop')
        op1 = Operator(eqn, dse='aggressive')

        op0(time_M=1)
        u0_norm = norm(u)

        u._data_with_inhalo[:] = 0.
        op1(time_M=1)
        u1_norm = norm(u)

        assert u0_norm == u1_norm

grid = Grid(shape=(8, 8))
        x, y = grid.dimensions
        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.

        eqn = Eq(u.forward, ((u[t, x, y] + u[t, x+1, y+1])*3*f +
                             (u[t, x+2, y+2] + u[t, x+3, y+3])*3*f + 1))
        op0 = Operator(eqn, dse='noop')
        op1 = Operator(eqn, dse='aggressive')

        op0(time_M=1)
        u0_norm = norm(u)

        u._data_with_inhalo[:] = 0.
        op1(time_M=1)
        u1_norm = norm(u)

        assert u0_norm == u1_norm

"""
    day = np.datetime64('today')
    try:
        l = np.load("norms%s.npy" % len(shape), allow_pickle=True)
        assert l[-1] == day
    except:
        tn = 500.  # Final time
        nrec = 130  # Number of receivers

        # Create solver from preset
        solver = acoustic_setup(shape=shape, spacing=[15. for _ in shape],
                                tn=tn, space_order=so, nrec=nrec,
                                preset='layers-isotropic', dtype=np.float64)
        # Run forward operator
        rec, u, _ = solver.forward()
        Eu = norm(u)
        Erec = norm(rec)

        # Run adjoint operator
        srca, v, _ = solver.adjoint(rec=rec)
        Ev = norm(v)
        Esrca = norm(srca)

        np.save("norms%s.npy" % len(shape), (Eu, Erec, Ev, Esrca, day), allow_pickle=True)

"""
        Computes the norms of the outputs in serial mode to compare with
        """
        shape = self._shapes[nd]
        so = self._so[nd]
        tn = 500.  # Final time
        nrec = 130  # Number of receivers

        # Create solver from preset
        solver = acoustic_setup(shape=shape, spacing=[15. for _ in shape],
                                tn=tn, space_order=so, nrec=nrec,
                                preset='layers-isotropic', dtype=np.float64)
        # Run forward operator
        rec, u, _ = solver.forward()
        Eu = norm(u)
        Erec = norm(rec)

        # Run adjoint operator
        srca, v, _ = solver.adjoint(rec=rec)
        Ev = norm(v)
        Esrca = norm(srca)

        return Eu, Erec, Ev, Esrca

m.data[:] = 1./(v*v)

        pde = m * u.dt2 - u.laplace
        stencil = Eq(u.forward, solve(pde, u.forward))

        src_term = src.inject(field=u.forward, expr=src * dt**2 / m)
        rec_term = rec.interpolate(expr=u.forward)

        op = Operator([stencil] + src_term + rec_term)

        # Make sure we've indeed generated OpenACC code
        assert 'acc parallel' in str(op)

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

        assert np.isclose(norm(rec), 490.56, atol=1e-2, rtol=0)

first_derivative(Gz * sintheta * sinphi, dim=y, fd_order=soh, matvec=T) +
               first_derivative(Gz * costheta, dim=z, fd_order=soh, matvec=T))

        # Equation
        eqn = [Eq(r.backward, Gzz)]

        op0 = Operator(eqn, subs=grid.spacing_map, opt=('noop', {'openmp': True}))
        op1 = Operator(eqn, subs=grid.spacing_map,
                       opt=('advanced', {'openmp': True, 'cire-mincost-sops': 1}))

        # Check numerical output
        op0(time_M=1)
        exp_norm_r = norm(r)
        r.data_with_halo[:] = 0.5
        summary = op1(time_M=1)
        assert np.isclose(norm(r), exp_norm_r, rtol=1e-5)

        # Also check against expected operation count to make sure
        # all redundancies have been detected correctly
        assert summary[('section1', None)].ops == 39

def test_viscoelastic():
    _, _, _, [rec1, rec2, v, tau] = run()
    assert np.isclose(norm(rec1), 12.28040, atol=1e-3, rtol=0)
    assert np.isclose(norm(rec2), 0.312461, atol=1e-3, rtol=0)

d_syn = Receiver(name='d_syn', grid=model.grid,
                         time_range=geometry.time_axis,
                         coordinates=geometry.rec_positions)
        # Update source location
        solver.geometry.src_positions[0, :] = source_locations[i, :]

        # Generate synthetic data from true model
        solver.forward(vp=model.vp, rec=d_obs)

        # Compute smooth data and full forward wavefield u0
        _, u0, _ = solver.forward(vp=vp_in, save=True, rec=d_syn)

        # Compute gradient from data residual and update objective function
        residual = compute_residual(residual, d_obs, d_syn)

        objective += .5*norm(residual)**2
        solver.gradient(rec=residual, u=u0, vp=vp_in, grad=grad)

    return objective, grad

source_mask_f = TimeFunction(name='source_mask_f', grid=model.grid,
                             time_order=0, dtype=np.int32)

source_mask_f.data[0, :, :, :] = source_mask.data[:, :, :]

eq1 = Eq(zind, source_id, implicit_dims=(time, x, y, z))
eq2 = Inc(u2.forward, source_mask * save_src[time, zind])


op2 = Operator([eq1, eq2])
op2.apply()
print(op2.ccode)

print(norm(u))
print(norm(u2))
import pdb; pdb.set_trace()
assert np.isclose(norm(u), norm(u2), atol=1e-06)

op2.apply()
print("Norm(u2):", norm(u2))
print("-----")


print("Norm(u):", norm(u))

print("Norm(u2):", norm(u2))

print("Norm(uref):", norm(uref))

print(norm(uref))

# import pdb; pdb.set_trace()

assert np.isclose(norm(uref), norm(u2), atol=1e-06)