Source:¶

Also thanks to Stephan Hoyer @shoyer for help with https://github.com/pydata/xarray/issues/1773

In [1]:
from sympy import *
init_printing()

import numpy as np
import xarray as xr
import itertools as Iter

import yt
yt.toggle_interactivity()

from scipy import constants as consts

import matplotlib.pyplot as plt
%matplotlib notebook

Using matplotlib backend: Qt5Agg
Warning: Cannot change to a different GUI toolkit: notebook. Using qt instead.


Sympy Aspect¶

In [2]:
I, x, y, z, r, rho0, theta=symbols('I, x, y, z, r, rho_0, theta')

In [3]:
rdef=sqrt(x**2+y**2+z**2)

Out[3]:
$$\left ( \sqrt{x^{2} + y^{2} + z^{2}}, \quad \operatorname{atan_{2}}{\left (z,\rho_{0} \right )}\right )$$
In [4]:
Ay=(consts.mu_0*rho0**2*I/4)*sin(theta)/r**2
Ay

Out[4]:
$$\frac{I \rho_{0}^{2}}{r^{2}} 3.14159265358979 \cdot 10^{-7} \sin{\left (\theta \right )}$$
In [5]:
Ay=simplify(Ay.subs({theta:thetaDef, r: rdef})); Ay

Out[5]:
$$\frac{3.14159265358979 \cdot 10^{-7} I \rho_{0}^{2} z}{\sqrt{\rho_{0}^{2} + z^{2}} \left(x^{2} + y^{2} + z^{2}\right)}$$
In [6]:
sups={rho0:.5, I:1} #current loop with radidis .5m with 1A
sups

Out[6]:
$$\left \{ I : 1, \quad \rho_{0} : 0.5\right \}$$
In [7]:
AyN=lambdify((x, y, z), Ay.subs(sups), dummify=False)
AyN(1,1,1)

Out[7]:
$$2.34160491035e-08$$
In [8]:
Bx=simplify(-Derivative(Ay, z).doit()); Bx

Out[8]:
$$\frac{I \rho_{0}^{2}}{\left(\rho_{0}^{2} + z^{2}\right)^{\frac{3}{2}} \left(x^{2} + y^{2} + z^{2}\right)^{2}} \left(6.28318530717959 \cdot 10^{-7} z^{2} \left(\rho_{0}^{2} + z^{2}\right) + 3.14159265358979 \cdot 10^{-7} z^{2} \left(x^{2} + y^{2} + z^{2}\right) - 3.14159265358979 \cdot 10^{-7} \left(\rho_{0}^{2} + z^{2}\right) \left(x^{2} + y^{2} + z^{2}\right)\right)$$
In [9]:
BxN=lambdify((x, y, z), Bx.subs(sups), dummify=False)
BxN(1,1,1)

Out[9]:
$$1.0927489581618925e-08$$
In [10]:
Bz=simplify(Derivative(Ay, x).doit()); Bz

Out[10]:
$$- \frac{6.28318530717959 \cdot 10^{-7} I \rho_{0}^{2} x z}{\sqrt{\rho_{0}^{2} + z^{2}} \left(x^{2} + y^{2} + z^{2}\right)^{2}}$$
In [11]:
BzN=lambdify((x, y, z), Bz.subs(sups), dummify=False)
BzN(1,1,1)

Out[11]:
$$-1.56106994023e-08$$

Generate the xarray dataset for the data¶

In [12]:
DomainSpaceSize = 6 # using cartesian 3D
SpaceDensity = [10] # 100 divisions in space 5 in time

In [13]:
x_coord = np.linspace(-DomainSpaceSize, +DomainSpaceSize, SpaceDensity[0])
y_coord = np.linspace(-DomainSpaceSize, +DomainSpaceSize, SpaceDensity[0])
z_coord = np.linspace(-DomainSpaceSize, +DomainSpaceSize, SpaceDensity[0])

In [14]:
CoorEnterFunc = lambda x, y, z: 1+0*x+0*y+0*z

dx=xr.DataArray(x_coord, dims='x')
dy=xr.DataArray(y_coord, dims='y')
dz=xr.DataArray(z_coord, dims='z')

In [15]:
SimDataSet = xr.Dataset({'Ay':(['x', 'y', 'z'], AyN(*np.meshgrid(x_coord, y_coord, z_coord)))},
coords={'x':x_coord, 'y':y_coord, 'z':z_coord})
SimDataSet.attrs['units']={'Ay':'V/s/m'}
SimDataSet

Out[15]:
<xarray.Dataset>
Dimensions:  (x: 10, y: 10, z: 10)
Coordinates:
* x        (x) float64 -6.0 -4.667 -3.333 -2.0 -0.6667 0.6667 2.0 3.333 ...
* y        (y) float64 -6.0 -4.667 -3.333 -2.0 -0.6667 0.6667 2.0 3.333 ...
* z        (z) float64 -6.0 -4.667 -3.333 -2.0 -0.6667 0.6667 2.0 3.333 ...
Data variables:
Ay       (x, y, z) float64 -7.247e-10 -8.327e-10 -9.345e-10 -1.003e-09 ...
Attributes:
units:    {'Ay': 'V/s/m'}
In [16]:
SimDataSet['Bx']=(['x', 'y', 'z',],BxN(*np.meshgrid(x_coord, y_coord, z_coord)))
SimDataSet.attrs['units']['Bx']='T'
SimDataSet

Out[16]:
<xarray.Dataset>
Dimensions:  (x: 10, y: 10, z: 10)
Coordinates:
* x        (x) float64 -6.0 -4.667 -3.333 -2.0 -0.6667 0.6667 2.0 3.333 ...
* y        (y) float64 -6.0 -4.667 -3.333 -2.0 -0.6667 0.6667 2.0 3.333 ...
* z        (z) float64 -6.0 -4.667 -3.333 -2.0 -0.6667 0.6667 2.0 3.333 ...
Data variables:
Ay       (x, y, z) float64 -7.247e-10 -8.327e-10 -9.345e-10 -1.003e-09 ...
Bx       (x, y, z) float64 7.969e-11 8.085e-11 6.879e-11 2.328e-11 ...
Attributes:
units:    {'Ay': 'V/s/m', 'Bx': 'T'}
In [17]:
SimDataSet['Bz']=(['x', 'y', 'z',],BzN(*np.meshgrid(x_coord, y_coord, z_coord)))
SimDataSet.attrs['units']['Bz']='T'
SimDataSet

Out[17]:
<xarray.Dataset>
Dimensions:  (x: 10, y: 10, z: 10)
Coordinates:
* x        (x) float64 -6.0 -4.667 -3.333 -2.0 -0.6667 0.6667 2.0 3.333 ...
* y        (y) float64 -6.0 -4.667 -3.333 -2.0 -0.6667 0.6667 2.0 3.333 ...
* z        (z) float64 -6.0 -4.667 -3.333 -2.0 -0.6667 0.6667 2.0 3.333 ...
Data variables:
Ay       (x, y, z) float64 -7.247e-10 -8.327e-10 -9.345e-10 -1.003e-09 ...
Bx       (x, y, z) float64 7.969e-11 8.085e-11 6.879e-11 2.328e-11 ...
Bz       (x, y, z) float64 -8.052e-11 -1.066e-10 -1.349e-10 -1.583e-10 ...
Attributes:
units:    {'Ay': 'V/s/m', 'Bx': 'T', 'Bz': 'T'}

Transfer the xarray to yt¶

In [18]:
bbox=[]
for i in ['x', 'y', 'z']:
bbox.append([float(SimDataSet.coords[i].min()), float(SimDataSet.coords[i].max())])
bbox=np.array(bbox)
bbox

Out[18]:
array([[-6.,  6.],
[-6.,  6.],
[-6.,  6.]])
In [19]:
YTDataObj= dict(MagPoty = (np.array(SimDataSet['Ay']), 'g/m**3'),
MagFieldx = (np.array(SimDataSet['Bx']), 'g/m**3'),
MagFieldz = (np.array(SimDataSet['Bz']), 'g/m**3')
)
#would not except unit in SimDataSet.attrs['units']['ChargeDen'] had to fake with
#'g/m**3'

In [20]:
YTDataObj = yt.load_uniform_grid(YTDataObj, tuple(dict(SimDataSet.dims).values()), length_unit="m", bbox=bbox, nprocs=64)

yt : [INFO     ] 2018-02-21 16:29:44,209 Parameters: current_time              = 0.0
yt : [INFO     ] 2018-02-21 16:29:44,211 Parameters: domain_dimensions         = [10 10 10]
yt : [INFO     ] 2018-02-21 16:29:44,212 Parameters: domain_left_edge          = [-6. -6. -6.]
yt : [INFO     ] 2018-02-21 16:29:44,212 Parameters: domain_right_edge         = [ 6.  6.  6.]
yt : [INFO     ] 2018-02-21 16:29:44,214 Parameters: cosmological_simulation   = 0.0


Trying to do transfer function but giving up

In [21]:
tf = yt.ColorTransferFunction((SimDataSet['Ay'].min(), SimDataSet['Ay'].max()), grey_opacity=False)
tf.map_to_colormap(SimDataSet['Ay'].min(), SimDataSet['Ay'].max(), scale=15.0, colormap='algae')


Plotting the Magnetic Potential¶

In [22]:
MagPotyScene = yt.create_scene(YTDataObj, 'MagPoty')
source2 = MagPotyScene[0]
MagPotyScene.camera.set_width(MagPotyScene.quan(7, 'm'))

In [23]:
MagPotyScene.show(sigma_clip=3)
# the transfer function should have a defualt to the max min of the data cube
#being renderd

yt : [INFO     ] 2018-02-21 16:29:45,197 Rendering scene (Can take a while).
yt : [INFO     ] 2018-02-21 16:29:45,200 Creating volume
/home/iridium/anaconda3/lib/python3.6/site-packages/yt/units/yt_array.py:1293: RuntimeWarning: invalid value encountered in log10
out_arr = func(np.asarray(inp), out=out, **kwargs)
yt : [INFO     ] 2018-02-21 16:29:45,315 Creating transfer function
yt : [INFO     ] 2018-02-21 16:29:45,316 Calculating data bounds. This may take a while.  Set the TranferFunctionHelper.bounds to avoid this.


Plotting the Bx and Bz fields (attempting superposition from within yt)¶

In [24]:
from yt.visualization.volume_rendering.api import Scene, VolumeSource

In [25]:
sc = Scene()
cam.resolution = [400, 400]
cam.position = YTDataObj.arr([9, 9, 9], 'm')
cam.switch_orientation()

In [26]:
BFx = VolumeSource(YTDataObj, field='MagFieldx')
BFx.use_ghost_zones = True
sc.show( sigma_clip=1)

yt : [INFO     ] 2018-02-21 16:29:49,025 Rendering scene (Can take a while).
yt : [INFO     ] 2018-02-21 16:29:49,027 Creating volume
/home/iridium/anaconda3/lib/python3.6/site-packages/yt/units/yt_array.py:1293: RuntimeWarning: invalid value encountered in log10
out_arr = func(np.asarray(inp), out=out, **kwargs)
yt : [INFO     ] 2018-02-21 16:29:49,692 Creating transfer function
yt : [INFO     ] 2018-02-21 16:29:49,693 Calculating data bounds. This may take a while.  Set the TranferFunctionHelper.bounds to avoid this.

In [27]:
# add rendering of x-velocity field
BFz = VolumeSource(YTDataObj, field='MagFieldz')
BFz.use_ghost_zones = True