Exact solution used in MES runs

We would like to MES the operation (in a cylindrical geometry)

$$ \nabla \cdot \left(\mathbf{u}_E\cdot\nabla \left[\frac{\nabla_\perp \phi}{B}n\right]\right) $$

As we have a homogenenous $B$-field, we have normalized it out, and remain with

$$ \nabla \cdot \left(\mathbf{u}_E\cdot\nabla\left[n\nabla_\perp \phi\right]\right) $$

In [1]:
%matplotlib notebook

from IPython.display import display

from sympy import init_printing
from sympy import S, Function, Derivative, Eq
from sympy import symbols, simplify, sympify, collect, expand

from boutdata.mms import x, y, z, t

import os, sys
# If we add to sys.path, then it must be an absolute path
common_dir = os.path.abspath("./../../common/")
# Sys path is a list of system paths
from CELMAPy.MES.mesGenerator import get_metric



In [2]:
metric = get_metric()
phi = Function('phi')(x,z)
n   = Function('n')  (x,z)

In [3]:
def DDX(f):
    return Derivative(f, x)

def DDZ(f):
    return Derivative(f, z)

Define the variables

In [4]:
# Initialization
the_vars = {}

Define functions

One can show that in cylindrical geometry

$$\mathbf{u}_E\cdot\nabla f = \frac{1}{J}(\partial_\theta \phi \partial_\rho f - \partial_\rho \phi \partial_\theta f)$$

Further on, we have that

$$ \partial_\rho \mathbf{e}^\rho = 0\\ \partial_\theta \mathbf{e}^\rho = \rho\mathbf{e}^\theta\\ \partial_\rho \mathbf{e}^\theta = -\frac{\mathbf{e}^\theta}{\rho}\\ \partial_\theta \mathbf{e}^\theta = -\frac{\mathbf{e}^\rho}{\rho} $$

and that

$$ n\nabla_\perp \phi = \mathbf{e}^\rho n \partial_\rho \phi + \mathbf{e}^\theta n \partial_\theta \phi $$

This means that one of the components can be written

$$ n\nabla_\perp \phi =\mathbf{e}^j n \partial_j \phi $$

such that

$$ \partial_i \left( n\nabla_\perp \phi \right) = \mathbf{e}^j\partial_i \left(n \partial_j \phi \right) + n \partial_j \phi \partial_i \left(\mathbf{e}^j\right) $$

This gives

\begin{align} \partial_\rho \left( \left[ n\nabla_\perp \phi \right]_\rho \right) &= \mathbf{e}^\rho \partial_\rho \left(n \partial_\rho \phi \right) + n \partial_\rho \phi \partial_\rho \left(\mathbf{e}^\rho\right) \\ &= \mathbf{e}^\rho \partial_\rho \left(n \partial_\rho \phi \right) \end{align}


\begin{align} \partial_\rho \left( \left[ n\nabla_\perp \phi \right]_\theta \right) &= \mathbf{e}^\theta \partial_\rho \left(n \partial_\theta \phi \right) + n \partial_\theta \phi \partial_\rho \left(\mathbf{e}^\theta\right) \\ &= \mathbf{e}^\theta \partial_\rho \left(n \partial_\theta \phi \right) - \frac{\mathbf{e}^\theta}{\rho}n \partial_\theta \phi \partial_\rho \end{align}

so that

\begin{align} \partial_\rho \left( n\nabla_\perp \phi \right) &= \mathbf{e}^\rho \partial_\rho \left(n \partial_\rho \phi \right) + \mathbf{e}^\theta \partial_\rho \left(n \partial_\theta \phi \right) - \frac{\mathbf{e}^\theta}{\rho}n \partial_\theta \phi \\ &= \mathbf{e}^\rho \left( \partial_\rho \left(n \partial_\rho \phi \right) \right) + \mathbf{e}^\theta \left( \partial_\rho \left(n \partial_\theta \phi \right) - \frac{1}{\rho}n \partial_\theta \phi \right) \end{align}

Which means that

\begin{align} \frac{1}{J}\partial_\theta \phi \partial_\rho \left( n\nabla_\perp \phi \right) = \mathbf{e}^\rho \frac{1}{J}\partial_\theta \phi \left( \partial_\rho \left(n \partial_\rho \phi \right) \right) + \mathbf{e}^\theta \frac{1}{J}\partial_\theta \phi \left( \partial_\rho \left(n \partial_\theta \phi \right) - \frac{1}{\rho}n \partial_\theta \phi \right) \end{align}

Further on, we have that

\begin{align} \partial_\theta \left( \left[ n\nabla_\perp \phi \right]_\rho \right) &= \mathbf{e}^\rho \partial_\theta \left(n \partial_\rho \phi \right) + n \partial_\rho \phi \partial_\theta \left(\mathbf{e}^\rho\right) \\ &= \mathbf{e}^\rho \partial_\theta \left(n \partial_\rho \phi \right) + \rho\mathbf{e}^\theta n \partial_\rho \phi \end{align}


\begin{align} \partial_\theta \left( \left[ n\nabla_\perp \phi \right]_\theta \right) &= \mathbf{e}^\theta \partial_\theta \left(n \partial_\theta \phi \right) + n \partial_\theta \phi \partial_\theta \left(\mathbf{e}^\theta\right) \\ &= \mathbf{e}^\theta \partial_\theta \left(n \partial_\theta \phi \right) - \frac{\mathbf{e}^\rho }{\rho} n \partial_\theta \phi \end{align}

so that

\begin{align} \partial_\theta \left( n\nabla_\perp \phi \right) &= \mathbf{e}^\rho \partial_\theta \left(n \partial_\rho \phi \right) + \mathbf{e}^\theta \rho n \partial_\rho \phi + \mathbf{e}^\theta \partial_\theta \left(n \partial_\theta \phi \right) - \frac{\mathbf{e}^\rho}{\rho} n \partial_\theta \phi \\ &= \mathbf{e}^\rho \left( \partial_\theta \left(n \partial_\rho \phi \right) - \frac{1}{\rho} n \partial_\theta \phi \right) + \mathbf{e}^\theta \left( \rho n \partial_\rho \phi + \partial_\theta \left(n \partial_\theta \phi \right) \right) \end{align}

Which means that

\begin{align} \frac{1}{J}\partial_\rho \phi \partial_\theta \left( n\nabla_\perp \phi \right) = \mathbf{e}^\rho \frac{1}{J}\partial_\rho \phi \left( \partial_\theta \left(n \partial_\rho \phi \right) - \frac{1}{\rho} n \partial_\theta \phi \right) + \mathbf{e}^\theta \frac{1}{J}\partial_\rho \phi \left( \rho n \partial_\rho \phi + \partial_\theta \left(n \partial_\theta \phi \right) \right) \end{align}

Collecting elements gives

\begin{align} \mathbf{u}_E\cdot\nabla \left( n\nabla_\perp \phi \right) =& \mathbf{e}^\rho \frac{1}{J} \left[ \partial_\theta \phi \partial_\rho \left(n \partial_\rho \phi \right) - \partial_\rho \phi \left( \partial_\theta \left(n \partial_\rho \phi \right) - \frac{1}{\rho} n \partial_\theta \phi \right) \right] \\& + \mathbf{e}^\theta \frac{1}{J} \left[ \partial_\theta \phi \left( \partial_\rho \left(n \partial_\theta \phi \right) - \frac{1}{\rho}n \partial_\theta \phi \right) - \partial_\rho \phi \left( \rho n \partial_\rho \phi + \partial_\theta \left(n \partial_\theta \phi \right) \right) \right] \end{align}

In [5]:
rho_element = (1/metric.J)*(
                                - (1/metric.J)*n*DDZ(phi)
display(Eq(symbols('v_rho'), rho_element.doit()))

In [6]:
theta_element = (1/metric.J)*(
                               DDX(n * DDZ(phi))
                            - (1/metric.J) * n * DDZ(phi)
                            + DDZ(n*DDZ(phi)) 

display(Eq(symbols('v_theta'), theta_element.doit()))

We have that

$$ \nabla \cdot \mathbf{A} = \frac{1}{J}\partial_i (J A^i) $$

And that

$$ A^i = g^{ij}A_j $$


$$ g^{ij} = 0 \iff i\neq j\\ g^{\rho\rho} =1 \\ g^{\theta\theta} = \frac{1}{\rho^2} $$

This means that

$$ \nabla \cdot \mathbf{A}_\perp = \frac{1}{J}\partial_i \left(Jg^{ik}A_k\right) = \frac{1}{\rho}\partial_\rho \left(\rho A_\rho\right) + \frac{1}{\rho}\partial_\theta \left(\frac{1}{\rho}A_\rho\right) $$

so that

$$ \nabla \cdot \left(\mathbf{u}_E\cdot\nabla\left[n\nabla_\perp \phi\right]\right) $$

can be written

In [7]:
S = (1/metric.J)*(
    + DDZ((1/metric.J)*theta_element)

display(Eq(symbols('S'), S))

In [8]:
S = simplify(S.doit())

display(Eq(symbols('S'), S))

Well, that is a mess, we should collect it another manner

In [9]:
strS = str(S)

# phi x derivatives
strS = strS.replace('Derivative(phi(x, z), x)', 'phi_x')
strS = strS.replace('Derivative(phi(x, z), x, x)', 'phi_xx')
strS = strS.replace('Derivative(phi(x, z), x, x, x)', 'phi_xxx')
# phi z derivatives
strS = strS.replace('Derivative(phi(x, z), z)', 'phi_z')
strS = strS.replace('Derivative(phi(x, z), z, z)', 'phi_zz')
strS = strS.replace('Derivative(phi(x, z), z, z, z)', 'phi_zzz')
# phi mixed derivatives
strS = strS.replace('Derivative(phi(x, z), x, z)', 'phi_xz')
strS = strS.replace('Derivative(phi(x, z), x, z, z)', 'phi_xzz')
strS = strS.replace('Derivative(phi(x, z), x, x, z)', 'phi_xxz')
# Non-derivatives
strS = strS.replace('phi(x, z)', 'phi')

# n x derivatives
strS = strS.replace('Derivative(n(x, z), x)', 'n_x')
strS = strS.replace('Derivative(n(x, z), x, x)', 'n_xx')
# n z derivatives
strS = strS.replace('Derivative(n(x, z), z)', 'n_z')
strS = strS.replace('Derivative(n(x, z), z, z)', 'n_zz')
# n mixed derivatives
strS = strS.replace('Derivative(n(x, z), x, z)', 'n_xz')
# Non-derivatives
strS = strS.replace('n(x, z)', 'n')

newS = sympify(strS)

In [10]:

Printing for comparison

In [11]:
display(Eq(symbols('S'), expand(S)))

In [12]:
display(Eq(symbols('S_new'), expand(newS)))

Port to BOUT++

In [13]:
newerS = collect(newS, symbols('n, phi_x, phi_z'), exact=True)

In [14]:

In [15]:
# Manual simplification (copied from manual_simplification.txt)
manSim = (sympify(\
'(  n*(-phi_x*(phi_xxz*x**3  + phi_xz*x**2  + phi_z*x   + phi_zzz*x) + phi_z*(phi_xx*x**2 + phi_xxx*x**3 + phi_xzz*x - 2*phi_zz)) - phi_x**2*(n_xz*x**3 + n_z*x**2) + phi_z**2*(n_xz*x - n_z) + phi_x*(  phi_z*(n_x*x**2 + n_xx*x**3- n_zz*x) - 2*n_z*(phi_xx*x**3 + phi_zz*x)) + 2*n_x*phi_z*(phi_xx*x**3 + phi_zz*x))/x**4'

In [16]:
display(expand(newerS)- expand(manSim))