Exact solution used in MES runs

We would like to MES the operation

$$ \nabla \cdot \mathbf{f}_\perp $$

Using cylindrical geometry.



In [1]:

    
%matplotlib notebook

from sympy import init_printing
from sympy import S
from sympy import sin, cos, tanh, exp, pi, sqrt

from boutdata.mms import x, y, z, t
from boutdata.mms import Delp2, DDX, DDY, DDZ

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
sys.path.append(common_dir)
from CELMAPy.MES import get_metric, make_plot, BOUT_print

init_printing()

Initialize



In [2]:

    
folder = '../twoGaussians/'
metric = get_metric()

Define the variables



In [3]:

    
# Initialization
the_vars = {}

Define manifactured solutions

Due to orthogonality we have that

$$ S = \nabla \cdot \mathbf{f}_\perp = \nabla_\perp \cdot \mathbf{f}_\perp = \frac{1}{J} \partial_i \left(Jf^i\right) = \frac{1}{J} \partial_x \left(Jf^x\right) + \frac{1}{J} \partial_z \left(Jf^z\right) $$

In cylindrical coordinates $J=\rho$, so this gives

$$ f = \frac{1}{\rho} \partial_\rho \left(\rho f^\rho\right) + \frac{1}{\rho} \partial_\theta \left(\rho f^\theta\right) = \partial_\rho f^\rho + \frac{f^\rho}{\rho} + \partial_\theta f^\theta $$

NOTE:

z must be periodic
The field $f(\rho, \theta)$ must be of class infinity in $z=0$ and $z=2\pi$
The field $f(\rho, \theta)$ must be single valued when $\rho\to0$
The field $f(\rho, \theta)$ must be continuous in the $\rho$ direction with $f(\rho, \theta + \pi)$
Eventual BC in $\rho$ must be satisfied



In [4]:

    
# We need Lx
from boututils.options import BOUTOptions
myOpts = BOUTOptions(folder)
Lx = eval(myOpts.geom['Lx'])



In [5]:

    
# Two gaussians
# NOTE: S actually looks good

# The skew sinus
# In cartesian coordinates we would like a sinus with with a wave-vector in the direction
# 45 degrees with respect to the first quadrant. This can be achieved with a wave vector
# k = [1/sqrt(2), 1/sqrt(2)]
# sin((1/sqrt(2))*(x + y))
# We would like 2 nodes, so we may write
# sin((1/sqrt(2))*(x + y)*(2*pi/(2*Lx)))
# Rewriting this to cylindrical coordinates, gives
# sin((1/sqrt(2))*(x*(cos(z)+sin(z)))*(2*pi/(2*Lx)))

# The gaussian
# In cartesian coordinates we would like
# f = exp(-(1/(2*w^2))*((x-x0)^2 + (y-y0)^2))
# In cylindrical coordinates, this translates to
# f = exp(-(1/(2*w^2))*(x^2 + y^2 + x0^2 + y0^2 - 2*(x*x0+y*y0) ))
#   = exp(-(1/(2*w^2))*(rho^2 + rho0^2 - 2*rho*rho0*(cos(theta)*cos(theta0)+sin(theta)*sin(theta0)) ))
#   = exp(-(1/(2*w^2))*(rho^2 + rho0^2 - 2*rho*rho0*(cos(theta - theta0)) ))

# A parabola
# In cartesian coordinates, we have
# ((x-x0)/Lx)^2
# Chosing this function to have a zero value at the edge yields in cylindrical coordinates
# ((x*cos(z)+Lx)/(2*Lx))^2


w = 0.8*Lx
rho0 = 0.3*Lx
theta0 = 5*pi/4
the_vars['f^x'] = sin((1/sqrt(2))*(x*(cos(z)+sin(z)))*(2*pi/(2*Lx)))*\
                exp(-(1/(2*w**2))*(x**2 + rho0**2 - 2*x*rho0*(cos(z - theta0)) ))*\
                ((x*cos(z)+Lx)/(2*Lx))**2
        
# The gaussian

# In cartesian coordinates we would like
# f = exp(-(1/(2*w^2))*((x-x0)^2 + (y-y0)^2))
# In cylindrical coordinates, this translates to
# f = exp(-(1/(2*w^2))*(x^2 + y^2 + x0^2 + y0^2 - 2*(x*x0+y*y0) ))
#   = exp(-(1/(2*w^2))*(rho^2 + rho0^2 - 2*rho*rho0*(cos(theta)*cos(theta0)+sin(theta)*sin(theta0)) ))
#   = exp(-(1/(2*w^2))*(rho^2 + rho0^2 - 2*rho*rho0*(cos(theta - theta0)) ))

w = 0.5*Lx
rho0 = 0.2*Lx
theta0 = pi
the_vars['f^z'] = exp(-(1/(2*w**2))*(x**2 + rho0**2 - 2*x*rho0*(cos(z - theta0)) ))

Calculate the solution



In [6]:

    
the_vars['S'] =   (1/metric.J)*DDX(metric.J*the_vars['f^x'], metric=metric)\
                + 0\
                + (1/metric.J)*DDZ(metric.J*the_vars['f^z'], metric=metric)

Plot



In [7]:

    
make_plot(folder=folder, the_vars=the_vars, plot2d=True, include_aux=False)

Print the variables in BOUT++ format



In [8]:

    
BOUT_print(the_vars, rational=False)