In [ ]:
import yt
dir = '/home/ychen/data/0only_0314_h1_nojiggle/'
ds = yt.load(dir+'MHD_Jet_nojiggle_hdf5_plt_cnt_0010')
In [ ]:
le = ds.domain_left_edge/16
re = ds.domain_right_edge/16
box = ds.box(le, re)
plot = yt.SlicePlot(ds, 'x', fields=['current_density_z'], data_source=box, width=((5, 'kpc'), (10, 'kpc')))
In [ ]:
plot.set_width(((3, 'kpc'), (6, 'kpc')))
plot.set_cmap('current_density_z', 'RdBu')
plot.set_zlim('current_density_z', -1E-14, 1E-14)
plot.set_log('current_density_z', True, linthresh=1E-17)
plot.annotate_grids()
plot.show()
In [ ]:
import numpy as np
from yt.fields.derived_field import ValidateSpatial
from yt.funcs import just_one
from yt.units import dimensions
c = yt.physical_constants.c
mu_0 = yt.physical_constants.mu_0
ftype = 'gas'
unit_system = ds.unit_system
sl_left = slice(None, -2, None)
sl_right = slice(2, None, None)
div_fac = 2.0
sl_center = slice(1, -1, None)
j_factors = {dimensions.magnetic_field_cgs/dimensions.length: 4.0*np.pi/c,
dimensions.magnetic_field_mks/dimensions.length: mu_0}
def _current_density_x(field, data):
f = (data[ftype, "magnetic_field_z"][sl_center,sl_right,sl_center] -
data[ftype, "magnetic_field_z"][sl_center,sl_left,sl_center]) \
/ (div_fac*just_one(data["index", "dy"]))
f -= (data[ftype, "magnetic_field_y"][sl_center,sl_center,sl_right] -
data[ftype, "magnetic_field_y"][sl_center,sl_center,sl_left]) \
/ (div_fac*just_one(data["index", "dz"]))
new_field = data.ds.arr(np.zeros_like(data[ftype, "magnetic_field_z"],
dtype=np.float64), f.units)
new_field[sl_center, sl_center, sl_center] = f
return new_field/j_factors[new_field.units.dimensions]
def _current_density_y(field, data):
f = (data[ftype, "magnetic_field_x"][sl_center,sl_center,sl_right] -
data[ftype, "magnetic_field_x"][sl_center,sl_center,sl_left]) \
/ (div_fac*just_one(data["index", "dz"]))
f -= (data[ftype, "magnetic_field_z"][sl_right,sl_center,sl_center] -
data[ftype, "magnetic_field_z"][sl_left,sl_center,sl_center]) \
/ (div_fac*just_one(data["index", "dx"]))
new_field = data.ds.arr(np.zeros_like(data[ftype, "magnetic_field_z"],
dtype=np.float64), f.units)
new_field[sl_center, sl_center, sl_center] = f
return new_field/j_factors[new_field.units.dimensions]
def _current_density_z(field, data):
f = (data[ftype, "magnetic_field_y"][sl_right,sl_center,sl_center] -
data[ftype, "magnetic_field_y"][sl_left,sl_center,sl_center]) \
/ (div_fac*just_one(data["index", "dx"]))
f -= (data[ftype, "magnetic_field_x"][sl_center,sl_right,sl_center] -
data[ftype, "magnetic_field_x"][sl_center,sl_left,sl_center]) \
/ (div_fac*just_one(data["index", "dy"]))
new_field = data.ds.arr(np.zeros_like(data[ftype, "magnetic_field_z"],
dtype=np.float64), f.units)
new_field[sl_center, sl_center, sl_center] = f
return new_field/j_factors[new_field.units.dimensions]
curl_validators = [ValidateSpatial(1,
[(ftype, "magnetic_field_x"),
(ftype, "magnetic_field_y"),
(ftype, "magnetic_field_z")])]
# Determine the correct unit for the current density
if dimensions.current_mks in unit_system.base_units:
current_density_unit = unit_system["current_mks"]/unit_system["length"]**2
else:
current_density_unit = unit_system["current_cgs"]/unit_system["length"]**2
for ax in 'xyz':
n = "current_density_%s" % ax
ds.add_field((ftype, n), sampling_type="cell",
function=eval("_%s" % n),
units=current_density_unit,
validators=curl_validators,
force_override=True)
In [ ]:
data = ds.index.grids[100]
f = (data[ftype, "magnetic_field_y"][sl_right,sl_center,sl_center] -
data[ftype, "magnetic_field_y"][sl_left,sl_center,sl_center]) \
/ (div_fac*just_one(data["index", "dx"]))
f -= (data[ftype, "magnetic_field_x"][sl_center,sl_right,sl_center] -
data[ftype, "magnetic_field_x"][sl_center,sl_left,sl_center]) \
/ (div_fac*just_one(data["index", "dy"]))
new_field = data.ds.arr(np.zeros_like(data[ftype, "magnetic_field_z"],
dtype=np.float64), f.units)
In [ ]:
j_factors
In [ ]:
le = ds.domain_left_edge/16
re = ds.domain_right_edge/16
box = ds.box(le, re)
plot = yt.SlicePlot(ds, 'x', fields=['current_density_z'], data_source=box)
In [ ]:
plot.zoom(32)
plot.show()
In [ ]:
plot = yt.SlicePlot(ds, 'x', fields=['current_density_x'], data_source=box).zoom(32)
plot.show()
In [ ]:
plot = yt.SlicePlot(ds, 'x', fields=['current_density_y'], data_source=box).zoom(32)
plot.show()
In [ ]:
plot = yt.SlicePlot(ds, 'x', fields=['current_density_magnitude'], data_source=box).zoom(32)
plot.set_zlim('current_density_magnitude', 1E-18, 1E-13)
plot.show()
In [ ]:
import yt
import numpy as np
from yt.frontends.stream.api import load_uniform_grid
ddims = (16,16,16)
data1 = {"magnetic_field_x":(np.random.random(size=ddims),"T"),
"magnetic_field_y":(np.random.random(size=ddims),"T"),
"magnetic_field_z":(np.random.random(size=ddims),"T")}
data2 = {}
for field in data1:
data2[field] = (data1[field][0]*1E4, "gauss")
ds1 = load_uniform_grid(data1, ddims, unit_system="cgs")
ds2 = load_uniform_grid(data1, ddims, unit_system="mks")
In [ ]:
dd1 = ds1.all_data()
j_cgs = dd1['current_density_x']*dd1['dy']*dd1['dz']
j_cgs.in_units('statA')
In [ ]:
dd2 = ds2.all_data()
j_mks = dd2['current_density_x']*dd2['dy']*dd2['dz']
j_mks.in_units('A')
In [ ]:
j_mks.to_equivalent('statA', 'CGS')
In [ ]:
dd1['x'].in_units('cm')
In [ ]:
dd2['x'].in_units('m')
In [ ]:
from yt.testing import assert_almost_equal
assert_almost_equal(j_cgs.to_equivalent('A', 'SI'), j_mks, 12)
In [ ]:
assert_almost_equal(j_mks.to_equivalent('statA', 'CGS'), j_cgs, 3)
In [ ]:
from yt.units import dimensions
from yt.utilities.physical_constants import mu_0, c
j_factors = {dimensions.magnetic_field_cgs/dimensions.length: c*4.0*np.pi,
dimensions.magnetic_field_mks/dimensions.length: 1.0/mu_0}
In [ ]:
I = yt.units.A
In [ ]:
I.to_equivalent('statA', 'CGS')
In [ ]: