Divergence Operator

In this tutorial, we will show you how build a divergence operator with ConvDo.

The divergence operator maps a vector field to a scalar field. The divergence operator for a vector field \(\mathbf{u}=[u,v]\) is defined as follows:

\[ \nabla \cdot \mathbf{u}=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y} \]

For a flow filed which satisfies the Navier-Stokes equation, the divergence of the velocity field should be zero, i.e., \(\nabla \cdot \mathbf{u}=0\). Kolmogorov flow is a classical turbulent flow; a snapshot of the velocity is illustrated as follows:

import torch
import matplotlib.pyplot as plt

flow= torch.load('../../binaries/kf_flow.pt')
channel_names=["u","v"]

figs, axs = plt.subplots(1, 2)
for i, ax in enumerate(axs):
    ax.imshow(flow[i].numpy(),origin='lower')
    ax.set_title(channel_names[i])
plt.show()

To calculate the divergence of this flow field, we can simply build a Divergence operator:

from ConvDO import *

class Divergence(FieldOperations):

    def __init__(self, 
                 domain_u:Domain,
                 domain_v:Domain,
                 order:int, 
                 device="cpu", 
                 dtype=torch.float32) -> None:
        super().__init__(order, device, dtype)
        self.velocity=VectorValue(ScalarField(domain=domain_u),ScalarField(domain=domain_v))

    def __call__(self, u:torch.Tensor,v:torch.Tensor):
        self.velocity.ux.register_value(u)
        self.velocity.uy.register_value(v)
        return (self.nabla@self.velocity).value

The domain_u and domain_v should be instances of the Domain class, which gives the boundary condition of the fields.

The order here refers to the order of the finite difference scheme used to calculate the derivate. For periodic boundary conditions, supported orders include 2,4,6,8. For un-periodic boundaries, only a second-order scheme is available.

The device and dtype should be inconsistent with the tensor used for the operation.

We use the VectorValue class to generate the velocity; the component of the vector is two ScalarField initialized with the velocity tensor and the corresponding domain. The VectorValue class also supports numerical numbers as its component, and we can perform some simple operations on the VectorValue:

a=VectorValue(1,2)
b=VectorValue(3,4)
c=a+b
d=a-b
e=5*a
f=a@b
print(c,d,e,f)

(4,6) (-2,-2) (5,10) 11

self.nabla is the nabla operator, \(\nabla\). Other supported operator includes the Laplacian operator self.nabla2 = \(\nabla^2\) and gradient operators self.grad_x = \(\partial/\partial x\) and self.grad_y = \(\partial/\partial y\). For instance, the Divergence operator can also be written as

class Divergence_new(FieldOperations):

    def __init__(self, 
                 domain_u:Domain,
                 domain_v:Domain,
                 order:int, 
                 device="cpu", 
                 dtype=torch.float32) -> None:
        super().__init__(order, device, dtype)
        self.u=ScalarField(domain=domain_u)
        self.v=ScalarField(domain=domain_v)

    def __call__(self, u:torch.Tensor,v:torch.Tensor):
        self.u.register_value(u)
        self.v.register_value(v)
        return (self.grad_x*self.u+self.grad_y*self.v).value

If you want to use these operators out of the FieldOperations class, you can also simply create them as follows:

nabla = ConvNabla(order=2)
nabla2 = ConvLaplacian(order=2)
grad_x = ConvGrad(order=2,direction="x")
grad_y = ConvGrad(order=2,direction="y")

Then, we can set up the computational domain of the field and calculate the divergence. The boundary condition of Kolmogorov flow is periodic thus we will set up a periodic domain. Here, delta_x and delta_y is the physical distance between mesh cells in the field. Besides PeriodicDomain, you can directly use the Domain class to add different boundary conditions and even obstacles inside the field.

Now, let's calculate the divergence!

u=flow[0]
v=flow[1]
domain=PeriodicDomain(delta_x=2*3.14/64,delta_y=2*3.14/64)
divergence=Divergence(order=6,domain_u=domain,domain_v=domain)
residual=divergence(u,v).numpy().squeeze().squeeze()
plt.imshow(residual,origin="lower")
plt.title("Divergence")
plt.colorbar()
plt.show()

Since the operation is based on the convolution provided by PyTorch, it supports batched input and output. If your input of a ScalarField only contains 2 dimensions \(H \times W\), it will automatically convert to a new tensor with the shape of \(B\times C \times H \times W\). The shape of the output is always \(B\times C \times H \times W\).