<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Calculation of the divergence of the advection of the perpendicular gradient of the potential times density using Clebsch coordinates
We would here like to calculate
$$
\nabla\cdot\left(\mathbf{u}E\cdot\nabla\left[n \nabla\perp\phi \right]\right)
$$
using cylindrical Clebsch coordinates, as this tensor identity has not been found in the literature.
NOTE
Step1:  Calculation of the $E\times B$ advection
We would now like to calculate
$$
\zeta = \nabla\cdot\left(\mathbf{u}E \cdot\nabla\left[n\nabla\perp\phi\right]\right)
$$
We will do this by
Calculate $n\nabla_\perp\phi$
Define $\mathbf{u}_E$
By first calculating $\nabla_\perp\phi$
Calculate $\mathbf{u}E\cdot\nabla \left(n\nabla\perp\phi\right)$
To check the different contributions we also
Calculate $\mathbf{u}_E\cdot\nabla f$
Calculate $\mathbf{a}\cdot \left(n\nabla_\perp\phi\right)$
Take the divergence of the resulting vector
Compare this with $B{\phi,\Omega^D}$
Calculation of $n\nabla_\perp\phi$
Step2:  Defining $\mathbf{u}_E$
We have that 
$${u}E = - \frac{\nabla\perp\phi\times\mathbf{b}}{B}$$
Remember that we are working with normalized equations, so $B$ (which in reality is $\tilde{B}$) is equal to $1$.
NOTE
Step3:  NOTE
Step4:  Calculation of $\mathbf{u}E\cdot\nabla \left(n\nabla\perp\phi\right)$
Calculation of $\mathbf{u}_E\cdot\nabla f$
Step5:  Calculation of  $\mathbf{a}\cdot\nabla \left(n\nabla_\perp\phi\right)$
Step6:  Using covariant vector
Step7:  Using contravariant vector
Step8:  Calculation of full  $\mathbf{u}E\cdot\nabla \left(n\nabla\perp\phi\right)$
Step9:  Calculation of $\nabla\cdot\left(\mathbf{u}E\cdot\nabla\left[n\nabla\perp\phi\right]\right)$
Step10:  Comparison with $B{\phi,\Omega^D}$
In cylindrical Clebsch coordinates, we have that $\mathbf{u}_E\cdot\nabla = {\phi,\cdot}$. However, we have normalized our equations so that $\tilde{B}=1$. As $B$ from the Clebsch system is not constant, we can achieve normalization by multiplying the Poisson bracket with the un-normalized $B$ (from the Clebsch system).
We define the vorticity-like field $\Omega^D$ to be $\Omega^D = \nabla\cdot\left(n\nabla_\perp\phi\right)$. In the Clebsch system this is written as
Step11:  We now write $\chi = B{\phi,\Omega^D}$
Step12:  The difference $\epsilon$ between $\zeta = \nabla\cdot\left(\mathbf{u}E\cdot\nabla\left[n\nabla\perp\phi\right]\right)$ and $\chi = B{\phi,\Omega^D}$ is given by
$$\epsilon = \zeta - \chi$$
Step13:  In fact we see that
\begin{align}
  \epsilon 
  - \left(
      -\frac{1}{\rho}[\partial_\rho\phi]{n, \partial_\rho\phi}
      -\frac{1}{\rho^3}[\partial_\theta\phi]{n, \partial_\theta\phi}
      +\frac{1}{\rho^4}[\partial_\theta n][\partial_\theta\phi]^2
    \right)
    =\
  \epsilon 
  - \left(
      \frac{1}{\rho}[\partial_\rho\phi]{\partial_\rho\phi,n}
      +\frac{1}{\rho^3}[\partial_\theta\phi]{\partial_\theta\phi, n}
      +\frac{1}{\rho^4}[\partial_\theta n][\partial_\theta\phi]^2
    \right)
    =  
\end{align}
Step14:  What is more interesting is in fact that
\begin{align}
  \epsilon 
  - \xi
    =
  \epsilon 
  - \frac{B}{2}{\mathbf{u}_E\cdot\mathbf{u}_E, n}
\end{align}
Step15:  Where
\begin{align}
  \mathbf{u}_E\cdot\mathbf{u}_E
  =
\end{align}
Step16:  Note that the last term $\frac{1}{\rho^4}(\partial_\theta n)(\partial_\theta\phi)^2$ does not appear to come from the Poisson bracket. This is however the case, and comes from the part which contains 
$\frac{1}{2}\partial_\rho\left(\frac{1}{\rho}\partial_\theta \phi\right)^2 =
\left(\frac{1}{\rho}\partial_\theta \phi\right)\partial_\rho\left(\frac{1}{\rho}\partial_\theta \phi\right)$
as
$\partial_i (fg) = f \partial_i g + g \partial_i f$
To summarize, we have
\begin{align}
    \zeta - (\chi + \xi) =
\end{align}
Step17:  Printing for comparison
<ASSISTANT_TASK:>
Python Code:
from IPython.display import display
from sympy import symbols, simplify, sympify, expand
from sympy import init_printing
from sympy import Eq, Function
from clebschVector import ClebschVec
from clebschVector import div, grad, gradPerp, advVec
from common import rho, theta, poisson
from common import displayVec
init_printing()
u_z = symbols('u_z', real = True)
# In reality this is a function, but as it serves only as a dummy it is here defined as a symbol
# This makes it easier to replace
f = symbols('f', real = True)
phi = Function('phi')(rho, theta)
n   = Function('n')(rho, theta)
# Symbols for printing
zeta, chi, epsilon = symbols('zeta, chi, epsilon')
Explanation: Calculation of the divergence of the advection of the perpendicular gradient of the potential times density using Clebsch coordinates
We would here like to calculate
$$
\nabla\cdot\left(\mathbf{u}E\cdot\nabla\left[n \nabla\perp\phi \right]\right)
$$
using cylindrical Clebsch coordinates, as this tensor identity has not been found in the literature.
NOTE: These are normalized equations. As $B$ is constant, we can choose $B_0$ so that the normalized $\tilde{B}=1$, thus, $B$ is excluded from these equations.
Also, we would like to compare this with
$$
B{\phi,\Omega^D}
$$
End of explanation
nGradPerpPhi = gradPerp(phi)*n
displayVec(nGradPerpPhi)
Explanation: Calculation of the $E\times B$ advection
We would now like to calculate
$$
\zeta = \nabla\cdot\left(\mathbf{u}E \cdot\nabla\left[n\nabla\perp\phi\right]\right)
$$
We will do this by
Calculate $n\nabla_\perp\phi$
Define $\mathbf{u}_E$
By first calculating $\nabla_\perp\phi$
Calculate $\mathbf{u}E\cdot\nabla \left(n\nabla\perp\phi\right)$
To check the different contributions we also
Calculate $\mathbf{u}_E\cdot\nabla f$
Calculate $\mathbf{a}\cdot \left(n\nabla_\perp\phi\right)$
Take the divergence of the resulting vector
Compare this with $B{\phi,\Omega^D}$
Calculation of $n\nabla_\perp\phi$
End of explanation
# The basis-vectors are contravariant => components are covariant
eTheta = ClebschVec(rho=0, theta=1, z=0, covariant=True)
eRho = ClebschVec(rho=1, theta=0, z=0, covariant=True)
B = eTheta^eRho
displayVec(B, 'B')
Blen = B.len()
display(Eq(symbols('B'), Blen))
b = B/(B.len())
displayVec(b, 'b')
Explanation: Defining $\mathbf{u}_E$
We have that 
$${u}E = - \frac{\nabla\perp\phi\times\mathbf{b}}{B}$$
Remember that we are working with normalized equations, so $B$ (which in reality is $\tilde{B}$) is equal to $1$.
NOTE: It migth appear that there is a discrepancy between having a coordinate system where $B$ is not constant where we have derived equation where $B$ is constant. This is because the cylindrical coordinate system is not a Clebsch system, but the metrics coinside. The Poisson bracket is the only place where $B$ is used explicitly, and care must be taken. The workaround is easy: Just multiply the Poisson bracket with $B$ to make it correct in cylindrical coordinates.
End of explanation
gradPerpPhi = gradPerp(phi)
displayVec(gradPerpPhi)
# Normalized B
BTilde = 1
# Defining u_E
ue = - ((gradPerpPhi^b)/BTilde)
displayVec(ue, 'u_E')
Explanation: NOTE: Basis vectors in $B$ are covariant, so components are contravariant
Calculation of $\nabla_\perp\phi$
End of explanation
ueDotGrad_f = ue*grad(f)
display(ueDotGrad_f)
Explanation: Calculation of $\mathbf{u}E\cdot\nabla \left(n\nabla\perp\phi\right)$
Calculation of $\mathbf{u}_E\cdot\nabla f$
End of explanation
aRho,  aZ,  aTheta  = symbols('a^rho, a^z, a^theta')
a_Rho, a_Z, a_Theta = symbols('a_rho, a_z, a_theta')
aCov = ClebschVec(rho = a_Rho, z=a_Z, theta = a_Theta, covariant=True)
aCon = ClebschVec(rho = aRho,  z=aZ,  theta = aTheta,  covariant=False)
Explanation: Calculation of  $\mathbf{a}\cdot\nabla \left(n\nabla_\perp\phi\right)$
End of explanation
aCovDotNablaGradPhi = advVec(aCov, nGradPerpPhi)
displayVec(aCovDotNablaGradPhi)
Explanation: Using covariant vector
End of explanation
aConDotNablaGradPhi = advVec(aCon, nGradPerpPhi)
displayVec(aConDotNablaGradPhi)
Explanation: Using contravariant vector
End of explanation
ueDotGradnGradPerpPhi = advVec(ue, nGradPerpPhi)
displayVec(ueDotGradnGradPerpPhi.doitVec())
displayVec(ueDotGradnGradPerpPhi.doitVec().simplifyVec())
Explanation: Calculation of full  $\mathbf{u}E\cdot\nabla \left(n\nabla\perp\phi\right)$
End of explanation
div_ueDotGradnGradPerpPhi = div(ueDotGradnGradPerpPhi)
zetaFunc = div_ueDotGradnGradPerpPhi.doit().expand()
display(Eq(zeta, simplify(zetaFunc)))
Explanation: Calculation of $\nabla\cdot\left(\mathbf{u}E\cdot\nabla\left[n\nabla\perp\phi\right]\right)$
End of explanation
vortD = div(gradPerp(phi)*n)
display(Eq(symbols('Omega^D'), vortD.doit().expand()))
Explanation: Comparison with $B{\phi,\Omega^D}$
In cylindrical Clebsch coordinates, we have that $\mathbf{u}_E\cdot\nabla = {\phi,\cdot}$. However, we have normalized our equations so that $\tilde{B}=1$. As $B$ from the Clebsch system is not constant, we can achieve normalization by multiplying the Poisson bracket with the un-normalized $B$ (from the Clebsch system).
We define the vorticity-like field $\Omega^D$ to be $\Omega^D = \nabla\cdot\left(n\nabla_\perp\phi\right)$. In the Clebsch system this is written as
End of explanation
poissonPhiVortD = Blen*poisson(phi, vortD)
chiFunc = poissonPhiVortD.doit().expand()
display(Eq(chi, chiFunc))
Explanation: We now write $\chi = B{\phi,\Omega^D}$
End of explanation
epsilonFunc = (zetaFunc - chiFunc).expand()
display(Eq(epsilon, epsilonFunc))
Explanation: The difference $\epsilon$ between $\zeta = \nabla\cdot\left(\mathbf{u}E\cdot\nabla\left[n\nabla\perp\phi\right]\right)$ and $\chi = B{\phi,\Omega^D}$ is given by
$$\epsilon = \zeta - \chi$$
End of explanation
epsMinusCorrection = epsilonFunc\
       -\
       (\
        (1/rho)*phi.diff(rho)*poisson(phi.diff(rho), n)\
       +(1/(rho)**3)*phi.diff(theta)*poisson(phi.diff(theta),n)\
       +(1/(rho)**4)*n.diff(theta)*(phi.diff(theta))**2
       )
        
display(epsMinusCorrection.simplify())
Explanation: In fact we see that
\begin{align}
  \epsilon 
  - \left(
      -\frac{1}{\rho}[\partial_\rho\phi]{n, \partial_\rho\phi}
      -\frac{1}{\rho^3}[\partial_\theta\phi]{n, \partial_\theta\phi}
      +\frac{1}{\rho^4}[\partial_\theta n][\partial_\theta\phi]^2
    \right)
    =\
  \epsilon 
  - \left(
      \frac{1}{\rho}[\partial_\rho\phi]{\partial_\rho\phi,n}
      +\frac{1}{\rho^3}[\partial_\theta\phi]{\partial_\theta\phi, n}
      +\frac{1}{\rho^4}[\partial_\theta n][\partial_\theta\phi]^2
    \right)
    =  
\end{align}
End of explanation
xi = (Blen/2)*poisson(ue*ue, n).doit()
epsMinusNewCorr = epsilonFunc - (Blen/2)*poisson(ue*ue, n).doit()
display(epsMinusNewCorr.simplify())
Explanation: What is more interesting is in fact that
\begin{align}
  \epsilon 
  - \xi
    =
  \epsilon 
  - \frac{B}{2}{\mathbf{u}_E\cdot\mathbf{u}_E, n}
\end{align}
End of explanation
display((ue*ue).doit())
Explanation: Where
\begin{align}
  \mathbf{u}_E\cdot\mathbf{u}_E
  =
\end{align}
End of explanation
display((zetaFunc - (chiFunc + xi)).simplify())
Explanation: Note that the last term $\frac{1}{\rho^4}(\partial_\theta n)(\partial_\theta\phi)^2$ does not appear to come from the Poisson bracket. This is however the case, and comes from the part which contains 
$\frac{1}{2}\partial_\rho\left(\frac{1}{\rho}\partial_\theta \phi\right)^2 =
\left(\frac{1}{\rho}\partial_\theta \phi\right)\partial_\rho\left(\frac{1}{\rho}\partial_\theta \phi\right)$
as
$\partial_i (fg) = f \partial_i g + g \partial_i f$
To summarize, we have
\begin{align}
    \zeta - (\chi + \xi) =
\end{align}
End of explanation
S = expand(zetaFunc)
strS = str(S)
# phi rho derivatives
strS = strS.replace('Derivative(phi(rho, theta), rho)', 'phi_x')
strS = strS.replace('Derivative(phi(rho, theta), rho, rho)', 'phi_xx')
strS = strS.replace('Derivative(phi(rho, theta), rho, rho, rho)', 'phi_xxx')
# phi theta derivatives
strS = strS.replace('Derivative(phi(rho, theta), theta)', 'phi_z')
strS = strS.replace('Derivative(phi(rho, theta), theta, theta)', 'phi_zz')
strS = strS.replace('Derivative(phi(rho, theta), theta, theta, theta)', 'phi_zzz')
# phi mixed derivatives
strS = strS.replace('Derivative(phi(rho, theta), rho, theta)', 'phi_xz')
strS = strS.replace('Derivative(phi(rho, theta), rho, theta, theta)', 'phi_xzz')
strS = strS.replace('Derivative(phi(rho, theta), rho, rho, theta)', 'phi_xxz')
# Non-derivatives
strS = strS.replace('phi(rho, theta)', 'phi')
# n rho derivatives
strS = strS.replace('Derivative(n(rho, theta), rho)', 'n_x')
strS = strS.replace('Derivative(n(rho, theta), rho, rho)', 'n_xx')
# n theta derivatives
strS = strS.replace('Derivative(n(rho, theta), theta)', 'n_z')
strS = strS.replace('Derivative(n(rho, theta), theta, theta)', 'n_zz')
# n mixed derivatives
strS = strS.replace('Derivative(n(rho, theta), rho, theta)', 'n_xz')
# Non-derivatives
strS = strS.replace('n(rho, theta)', 'n')
newS = sympify(strS)
display(Eq(symbols('S_new'), expand(newS)))
Explanation: Printing for comparison
End of explanation
<END_TASK>