Simple second order accurate gradient

[cdelavegamartin]

Kind of related to this discussion , I've been implementing a second order accurate gradient, using centered differences in the interior, and one-sided second order scheme at the edges, ie, something where the matrix blocks for each component look like this for h=1:

[[-1.5  2.  -0.5  0.   0. ]
 [-0.5  0.   0.5  0.   0. ]
 [ 0.  -0.5  0.   0.5  0. ]
 [ 0.   0.  -0.5  0.   0.5]
 [ 0.   0.   0.5 -2.   1.5]]

With devito:

import devito as dv
import numpy as np

# grid
x_lims = (0,6)
y_lims = (0,4)
nx = 7
ny = 5
x = np.linspace(x_lims[0], x_lims[1], nx)
y = np.linspace(y_lims[0], y_lims[1], ny)
dx = x[1]-x[0]
dy = y[1]-y[0]
X,Y = np.meshgrid(x,y, indexing='ij')
u = np.sqrt((X)**2 + (Y)**2)
grid = dv.Grid(extent=(x_lims[1]-x_lims[0], y_lims[1]-y_lims[0]), shape=(nx, ny))
x,y = grid.dimensions

u_dv = dv.Function(name='u', grid=grid, space_order=3,)
# Initialize u_dv with the same data as u
u_dv.data[:] = u[:]
n_dv = dv.VectorFunction(name='n', grid=grid,space_order=2,staggered=(None,None))

#  Subdomains
class Middlex(dv.SubDomain):
    name = 'midx'
    def define(self, dimensions):
        x, y = dimensions
        return {x: ('middle', 1, 1), y: y}
    
class Middley(dv.SubDomain):
    name = 'midy'
    def define(self, dimensions):
        x, y = dimensions
        return {x: x, y: ('middle', 1, 1)}

class Left(dv.SubDomain):
    name = 'left'
    def define(self, dimensions):
        x, y = dimensions
        return {x: ('left', 1), y: y}
class Right(dv.SubDomain):
    name = 'right'
    def define(self, dimensions):
        x, y = dimensions
        return {x: ('right', 1), y: y}
class Bottom(dv.SubDomain):
    name = 'bottom'
    def define(self, dimensions):
        x, y = dimensions
        return {x: x, y: ('left', 1)}
class Top(dv.SubDomain):
    name = 'top'
    def define(self, dimensions):
        x, y = dimensions
        return {x: x, y: ('right', 1)}
midx = Middlex(grid=grid)
midy = Middley(grid=grid)
left = Left(grid=grid)
right = Right(grid=grid)
bottom = Bottom(grid=grid)
top = Top(grid=grid)

grad_u = dv.grad(u_dv, fd_order=3)
eq_nx = dv.Eq(n_dv[0], grad_u[0], subdomain=midx)
eq_ny = dv.Eq(n_dv[1], grad_u[1], subdomain=midy)
edge_eqs = []
edge_eqs.append(dv.Eq(n_dv[0], u_dv.dxr(fd_order=(2),),subdomain=left))
edge_eqs.append(dv.Eq(n_dv[0], u_dv.dxl(fd_order=(2),),subdomain=right))
edge_eqs.append(dv.Eq(n_dv[1], u_dv.dyr(fd_order=(2),),subdomain=bottom))
edge_eqs.append(dv.Eq(n_dv[1], u_dv.dyl(fd_order=(2),),subdomain=top))


eqs = [eq_nx,eq_ny] + edge_eqs
op = dv.Operator(eqs)

# Naive loop
n_loop = np.zeros((2,nx,ny))
for i in range(nx):
    for j in range(ny):
        if i==0:
            n_loop[0,i,j] = (-1.5*u[i,j]+2*u[i+1,j]-0.5*u[i+2,j])/dx
        elif i==nx-1:
            n_loop[0,i,j] = (1.5*u[i,j]-2*u[i-1,j]+0.5*u[i-2,j])/dx
        else:
            n_loop[0,i,j] = (u[i+1,j]-u[i-1,j])/(2*dx)
        
        if j==0:
            n_loop[1,i,j] = (-1.5*u[i,j]+2*u[i,j+1]-0.5*u[i,j+2])/dy
        elif j==ny-1:
            n_loop[1,i,j] = (1.5*u[i,j]-2*u[i,j-1]+0.5*u[i,j-2])/dy
        else:
            n_loop[1,i,j] = (u[i,j+1]-u[i,j-1])/(2*dy)

# compare
print(np.abs(n_dv[0].data - n_loop[0]).max())
print(np.abs(n_dv[1].data - n_loop[1]).max())

5.251605639244872e-07
5.892261540907384e-07

Questions:

• Is ther a more compact way to express this? specially when working in higher dimensions. I have also tried implementing the edges similar to the examples of boundary conditions, but it throws an error: TypeError: 'Derivative' object is not subscriptable

edge_eqs.append(dv.Eq(n_dv[0][0,y], u_dv.dxr(fd_order=(2),)[0,y]))
edge_eqs.append(dv.Eq(n_dv[0][nx-1,y], u_dv.dxl(fd_order=(2),)[0,y]))
edge_eqs.append(dv.Eq(n_dv[1][x,0], u_dv.dyr(fd_order=(2),)[0,y]))
edge_eqs.append(dv.Eq(n_dv[1][x,ny-1], u_dv.dyl(fd_order=(2),)[0,y]))

Or computes something that is incorrect

edge_eqs.append(dv.Eq(n_dv[0][0,y], u_dv.dxr(fd_order=(2),)))
edge_eqs.append(dv.Eq(n_dv[0][nx-1,y], u_dv.dxl(fd_order=(2),)))
edge_eqs.append(dv.Eq(n_dv[1][x,0], u_dv.dyr(fd_order=(2),)))
edge_eqs.append(dv.Eq(n_dv[1][x,ny-1], u_dv.dyl(fd_order=(2),)))

• what is the correct way to set the order of the gradient? setting order=2 doesn't work as expected, but fd_order=2 does ( Iassume is passing it to the child objects?

#  gradient fd_order
display(dv.grad(u_dv, fd_order=2).evaluate)
display(dv.grad(u_dv, fd_order=3).evaluate)
display(dv.grad(u_dv, order=2).evaluate)
display(dv.grad(u_dv, order=3).evaluate)