How to solve Time-dependent PDE with a dynamic boundary conditions in 2D disk?

[chunyanlimath]

Dear Prof.Lu and Communities,
I am trying to solve a time-dependent PDE equipped with a dynamic boundary condition in the 2D disk. Does anyone know how to implement such ad dynamic boundary condition?
The boundary condition is shown in the picture, which is a surface partial differential equation of u(x, t), where x \in \Omega, t \in[0, T], Ms, \beta, \alpha, are constant parameters. n is the unit outer normal vector of the boundary \partial \Omega. \nabla_s=\nabla - n(n\dot \nabla) is surface gradient operator.

As we can see that, there are several terms n dot nabla, and n dot nabla^3, n dot nabla_s and n dot nabla_s^3, which can't not be written as a non-homogeneous Neumann BC. Hence, I tried to define this boundary condition using OperatorBC, but I always get error
RuntimeError: Can't call numpy() on Tensor that requires grad. Use tensor.detach().numpy() instead. The first error shows in def func_boundary(x, u, _): normal = geomtime.boundary_normal(x[:, 0:2]) .

Here are the several ideas I tried but all failed.

1. I tried to replace this line by normal=geom.boundary_normal(x[:, 0:2]) resulted in the same error.
2. I tried to remove the time component, so the input is x =[x, y] in R^2, and output is u(x). And the time derivative in the boundary condition is also removed. But I still get the same error.
3. I tried to rewrite this boundary condition into an inhomogeneous Neumann type B.C. , but since there are several n dot nabla, n dot nanla^3, so it can't be written into that form.

Why this error comes out? Is it because I didn't correctly compute the normal vector of the boundary of a GeometryXTime Domain?
Any suggestions or hints will be appreciated!

Here is the code I wrote

# --------------- define computational geometry -----------
geom = dde.geometry.Disk([0, 0], R)`
timedomain = dde.geometry.TimeDomain(0, end_time)
geomtime = dde.geometry.GeometryXTime(geom, timedomain)
 
 # -------- define nabla u in bulk ----------- 
 
def nabla(x, u):   # x=[x, y, t] input, u is output of NN
    du_x = dde.grad.jacobian(u, x, i=0, j=0)
    du_y = dde.grad.jacobian(u, x, i=0, j=1)
    grad = np.array(
        [
            du_x,
            du_y
        ]
    )
    return grad   # shape=(2,)

 

# ------- define laplace u / 2nd derivative -----------
def nabla2(x, u):  # x=[x, y, t]
    du_xx = dde.grad.hessian(u, x, i=0, j=0)
    du_yy = dde.grad.hessian(u, x, i=1, j=1)
    return du_xx + du_yy   # scalar


# -------- define nabla3 u / 3rd derivative -----------
def nabla3(x, u):
    dnablau_x = dde.grad.jacobian(nabla2(x, u), x, i=0, j=0)
    dnablau_y = dde.grad.jacobian(nabla2(x, u), x, i=0, j=1)
    return np.array([dnablau_x, dnablau_y])   # shape=(2, ) # need check !!!


# ----------- define Laplace^2 u -------------
def nabla4(x, u):
    du_xx = dde.grad.hessian(u, x, i=0, j=0)  # H[i][j] = d^2y / dx_i dx_j, where i,j=0,…,dim_x-1.
    du_yy = dde.grad.hessian(u, x, i=1, j=1)
    du_xxxx = dde.grad.hessian(du_xx, x, i=0, j=0)
    du_yyyy = dde.grad.hessian(du_yy, x, i=1, j=1)
    du_xxyy = dde.grad.hessian(du_yy, x, i=0, j=0)
    return du_xxxx + du_yyyy + 2 * du_xxyy


# ----------- define mu in bulk ----------
def mu(x, u):
    return (1-eps) * u + u**3 + 2 * nabla2(x, u) + nabla4(x, u)


# --------- define n dot nabla --------------
def ndotnabla(x, u):
   # xx = x.detach().numpy()
    normal = geomtime.boundary_normal(x[:,0:2])  # only take the first 2 spatial components, the third one is time 
    normal = np.array([normal]).T
    return np.sum(nabla(x, u) * normal, axis=0)  # sum along row, keep column dim


# -------- define nabla_s u on surface -----------
def nabla_s(x, u):
    # xx = x.detach().numpy()
    normal = geomtime.boundary_normal(x[:, 0:2])
    normal = np.array([normal]).T
    return nabla(x, u) - normal * ndotnabla(x, u)


# -------- define nabla_s ^2 u on surface ------------
def nabla_2s(x, u):
    return nabla2(x, u) - np.sum(np.square(ndotnabla(x, u)), axis=0)


# ------ define nabla_s^3 u on surface -----------
def nabla_3s(x, u):
    return nabla_s(x,  nabla_2s(x, u))


# -------- define nabla_s ^4 u on surface ---------
def nabla_4s(x, u):
    return np.sum(nabla_s(x, nabla_3s(x, u)), axis=0)


# --------- define mu_s on surface -----------
def mu_s(x, u):
    # xx = x.detach().numpy()
    normal = geomtime.boundary_normal(x[:, 0:2])
    normal = np.array([normal]).T
    return (
            (1-eps)*u + u**3 - 2 * ndotnabla(x, u)
           - np.sum(normal * nabla_3s(x, u), axis=0) + 2*nabla_2s(x, u) + nabla_4s(x, u)
           + 4*H * np.sum(normal * nabla_s(x, u), axis=0)
           + 4*H**2 * nabla_2s(x, u) + 4 * H * np.sum(normal * nabla_3s(x, u), axis=0)
    )

# ------------ define boundary -------------
def bounary(x, on_boundary):
    return on_boundary    # and geomtime.on_boundary(x)

# ----------- define boundary value / eq ----------
def func_boundary(x, u, _):
    # xx = x.detach().numpy()
    normal = geomtime.boundary_normal(x[:, 0:2])
    normal = np.array([normal]).T
    du_t = dde.grad.jacobian(u, x, i=0, j=2)
    f = Ms2 * nabla_2s(x, mu_s(x, u)) - beta1**2 / alpha1 * mu_s(x, u) + beta1 / alpha1 * mu(x, u)
    return  du_t - f
# -------------- define B.C. ----------------
bc = dde.icbc.OperatorBC(geomtime, func_boundary, bounary)

[lululxvi]

You need to use tensorflow, not numpy.

[chunyanlimath]

Thanks for your comment, but my backend is pytorch. The way that I fixed it is to change the type of x (a tensor) to be a numpy array first, and then send it to geomtime.boundary_normal()function.