I am trying to write an example computation with multivariable total variation to include in my functional analysis notes using the following definition from Wikipedia:
Let $\Omega$ be an open subset of $\mathbb{R}^n$. Given a function $f$ belonging to $L^1(\Omega)$, the total variation of $f$ in $\Omega$ is defined as $$V(f,\Omega):=\sup \left \{ \int_\Omega f(x) \text{div} \phi(x) dx : \phi \in C_c^1(\Omega, \mathbb{R}^n), \|\phi\|_{L^\infty(\Omega)} \leq 1\right\}$$ $C_c^1(\Omega, \mathbb{R}^n)$ is the set of continuously differentiable vector functions of compact support contained in $\Omega$, $\|\cdot\|_{L^\infty(\Omega)}$ is the essential supremum norm, and $\text{div}$ is the divergence operator.
I would like to use this formula directly and demonstrate the process of taking the supremum.
I understand that if $f$ is $C^1$ on $\overline \Omega$, then the formula for total variation is simplified to the computation of
$$
V(f, \Omega) = \int_\Omega |\nabla f(x)| dx,
$$ which I am not trying to use here.
So, the problem is computing the total variation of $f: \mathbb{R}^2 \to \mathbb{R}$,
$$
f(x,y) =\frac{xy}{x^2 + y^2}
$$
on $\Omega$, where $\Omega$ is the open unit disk in $\mathbb{R}^2$. So, $\Omega = \{x : x \in \mathbb{R}^2 \text{ and }\|x\| < 1\}.$ Here is a visual of this situation:

My Attempted Solution.
$f$ is a classic example of a function discontinuous at 0, so $f \notin C^1(\overline \Omega)$. We first show that $f \in L^1(\Omega)$. Recall $f \in L^1(\Omega) \iff \int_\Omega |f| < \infty$. So, $$\begin{align} \int_{\Omega} |f| &= \iint_D |f(x,y)| dA \\ &= \int_{-1}^1 \int_{-\sqrt{1 - x^2}}^{\sqrt{1 - x^2}} \left \lvert \frac{xy}{x^2 + y^2}\right \rvert dy dx \\ &= 1 < \infty \qquad \text{(C.A.S)}\end{align}.$$ Hence, $f \in L^1(\Omega)$. It remains to compute $$V(f,\Omega):=\sup \left \{ \int_\Omega f(x,y) \text{div} \phi(x,y) dx : \phi \in C_c^1(\Omega, \mathbb{R}^2), \|\phi\|_{L^\infty(\Omega)} \leq 1\right\}$$
Before taking the supremum over $\phi \in C_c^1(\Omega, \mathbb{R}^2)$, we attempt the following simplification, $$\begin{align*} \int_\Omega f(x,y) \text{div} \phi(x,y) dx &= \iint_D f(x,y) \left(\nabla \cdot \left(\phi_x, \phi_y\right)\right) dA \\ &= \iint_D f(x,y) \left(\frac{\partial \phi_x}{\partial x} + \frac{\partial \phi_y}{\partial y}\right) dA \\ &= \iint_D f(x,y) \frac{\partial \phi_x}{\partial x} dA + \iint_D f(x,y) \frac{\partial \phi_y}{\partial y} dA \end{align*} $$
This is where I am stuck in terms of working with $\phi(x,y)$. It boils down to two main questions:
Given what we know about $\phi$, can the above be simplified any further to an expression that does not depend on $\phi$? (So that we don't have to take a supremum) Otherwise,
If we do have to take a supremum over the vector fields $\phi$, how would one go about this? I understand in the single variable case of total variation, one is simply constructing a family of partitions on an interval that enable you to take the supremum over all partitions. But here, how would I go about constructing a family of vector fields that enable me to take the supremum required?