Homec4science
Diffusion toma (master) Browse

R9690/postproc/spec_mag.py
f456046b0f48master

Last Change
Raw File

spec_mag.py
View Options

# Author : TOC
#
# generates plots of magnetic energy spectra for each pair of runs
# import modules
import pencil as pc
import numpy as np
import pylab as plt
import rundirs as rd
import fitfunc as ff
import os
from cycler import cycler
# latex fonts:
from matplotlib import rcParams
rcParams['text.usetex'] = True
rcParams['text.latex.preamble'] = [r'\usepackage{amssymb}']
# ``constants''
rhom = 1.
#Cmu = 16.
#Clambda = 1.
# cycle through run directories
for idx,val in enumerate(rd.runs):
os.chdir(val)
# read in data
par = pc.read.param(datadir='.', quiet=True)
p = pc.read.power(datadir='.')
eta = par.eta
#teta = 1/eta
lam5 = par.lambda5
mu50 = par.mu5_const
eos = par.lrelativistic_eos # True if relativistic
# calculate k_lambda from run parameters
#kl = np.sqrt(rhom*lam5*Cmu/Clambda) * eta * mu50
# print('k_lambda =', kl)
# this are the power spectra at different times
tspec = p.t[0:-1]
pmag = p.mag
minspec = 0
maxspec = len(tspec)
# print('length of tspec =', maxspec)
mintime = tspec[0]*eta
# print('mintime =', mintime)
maxtime = tspec[maxspec-1]*eta
# print('maxtime =', maxtime)
# karray = np.linspace(0, len(pmag[0]), num=len(pmag[0]))
plotspectra = 14
# find indicies for lin-distance in time:
# ideallinspace = np.linspace(tspec[0], tspec[maxspec-1],
# num=plotspectra)
# print('ideallinspace/teta =', np.array(ideallinspace)/teta)
# plotarray = []
# plotarrayindex = []
# j = 0
# jspec = 0
# for j in range(plotspectra):
# jspec = np.abs(tspec - ideallinspace[j]).argmin()
# plotarrayindex.append(jspec)
# plotarray.append(tspec[jspec])
# print('plotarrayindex =', plotarrayindex)
# print('Spectra times:', np.array(plotarray)/teta)
# to ensure corresponding nonrel and relativ runs on same plot
if idx%2==0:
# plot
fig = plt.figure()
ax = plt.subplot()
plt.clf()
plt.axis([1, 140, 1e-16, 1e0])
plt.title('Magnetic Energy Spectra for '+rd.names[int(idx/2)]+' Runs')
plt.ylabel('$E_\mathrm{mag}(k)$')
plt.xlabel('$k/k_1$')
plt.xscale('log')
plt.yscale('log')
# vertical dotted line at instability scale
plt.axvline(x=mu50/2, linestyle=':', lw=0.9)
# vertical dotted line at k_lambda (if exists)
#if 1<kl<20: plt.axvline(x=kl, linestyle=(':'), lw=0.9)
# k power lines
klin = np.linspace(2, 7, 1000)
plt.plot(klin, ff.powlaw(klin, a=1e-14, b=4), lw=0.5) # q=-4
plt.text(4, 1e-12, r'$k^4$', fontsize='small')
klin2 = np.linspace(8, 50, 2000)
plt.plot(klin2, ff.powlaw(klin2, a=10**1.5, b=-2), lw=0.5) # q=2
plt.text(22, 1e-1, r'$k^{-2}$', fontsize='small')
# klin3 = np.linspace(20, 70, 1000)
# klin4 = np.linspace(70, 300, 1000)
# plt.plot(klin3, ff.powlaw(klin3, a=20, b=-3), lw=0.5) # q=3
# plt.plot(klin4, ff.powlaw(klin4, a=1400, b=-4), lw=0.5) # q=4
# color bar (subplot 2)
a = np.array([[mintime, rd.spectimes[-1]]])
cmap = plt.cm.get_cmap('copper')
i = ax.imshow(a, cmap=cmap)
fig.colorbar(i, label='$t/t_\eta$')
# note that colorbar is a method of the figure, not the axes
n=0
# plot spectra (subplot 1)
# create line colors
colors = [plt.cm.copper(c) for c in np.linspace(0, 1, plotspectra)]
ax.set_prop_cycle(cycler('color', colors))
for i,v in enumerate(rd.spectimes): # plotarrayindex:
pspec = np.abs(tspec*eta - v).argmin()
#print('k max =', len(pmag[pspec]))
plt.plot(pmag[pspec], color=colors[n],
linewidth=0.8, linestyle=('--' if eos==True else '-'))
print('len(pmag[pspec]) =', len(pmag[pspec]))
n+=1
#print(min(pmag[pspec]),max(pmag[pspec]))
# spectra and horiz line at start of inverse cascade
#ichl = Cmu * rhom * mu50 * eta**2
for ind,vle in enumerate(pmag):
if np.argmax(vle)==19: # print('Em(kp) at IC:',max(vle))
#plt.axhline(y=max(vle), lw=0.8, color='g',
# linestyle=('-.' if eos==True else ':'))
plt.plot(vle, lw=0.8, color='g',
linestyle=('--' if eos==True else '-'))
#plt.axhline(y=ichl, lw=0.8, color='g',
# linestyle=('-.' if eos==True else ':'))
print('Cmu = Emic / (rhom*mu50*eta^2) =',
np.max(vle)/(rhom*mu50*(eta**2)))
break
#plt.axhline(y=ichl, lw=0.8, color='g',
# linestyle=('-.' if eos==True else ':'))
# spectra and horiz line at k_p = k_lambda :
#klhl = Clambda * mu50 / lam5
#if 1<=kl<=20:
# for ind1,vle1 in enumerate(pmag):
# if np.argmax(vle1)<=int(round(kl)):
#plt.axhline(y=max(vle1), lw=0.8, color='m',
# linestyle=('-.' if eos==True else ':'))
# plt.plot(vle1, lw=0.8, color='m',
# linestyle=('--' if eos==True else '-'))
#plt.axhline(y=klhl, lw=0.8, color='m',
# linestyle=('-.' if eos==True else ':'))
#print('t_lambda =', tspec[ind1]/teta)
# break
#plt.axhline(y=klhl, lw=0.8, color='m',
# linestyle=('-.' if eos==True else ':'))
# first spectra with k_p = 1 (box size)
Emmax = []
kpmax = []
Clk1 = 0
for ind2,vle2 in enumerate(pmag):
Emmax.append(np.max(vle2))
kpmax.append(np.argmax(vle2))
if np.argmax(vle2)==1:
plt.plot(vle2, color=colors[0], lw=0.8,
linestyle=('--' if eos==True else '-'))
#print('Clambda(kp=1) = Em * lambda5 / mu50 =',
Clk1 = np.max(vle2)*lam5/mu50 #)
break
ClEm = np.max(Emmax)*lam5/mu50 #, 'at k =', kpmax[np.argmax(Emmax)])
print('Clambda(Emmax) =', ClEm)
print('Clambda = Emsat * lambda5 / mu50 =', (Clk1+ClEm)/2,
'+/-', np.abs(ClEm - Clk1)/2)
for ind3,val3 in enumerate(tspec):
if val3*eta == rd.tsat[idx]:
print()
break
# to ensure corresponding nonrel and relativ runs on same plot
if idx%2==1:
os.chdir('/Users/TCB/Desktop/Research/JennySchober'
'/pencil-code/toma/postproc')
# save figure
fig.savefig('spec_mag_'+rd.els[int(idx/2)]+'.pdf',
bbox_inches='tight')
c4science · Help