diag_python/simple_example.py at main · GKV-developers/diag_python · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
#!/usr/bin/env python
# coding: utf-8

# ## Simple example of diag_python
# **In the below, lines from `import sys` to `geom_set` are initialization for data analysis.
# Users import a function of which he would like to analyze (`phiinxy` in this example).**

# In[1]:


### Append diag_python packages to path ###
import sys
sys.path.append("./src/")

### Read Zarr format data *.zarr/ by xarray ###
from diag_rb import rb_open
xr_phi = rb_open('../phi/gkvp.phi.*.zarr')

### Set geometric constants from GKV namelist etc. ###
from diag_geom import geom_set
geom_set(headpath='../src/gkvp_header.f90',
         nmlpath="../gkvp_namelist.001",
         mtrpath='../hst/gkvp.mtr.001')

# Plot phi[y,x] at t[it], zz[iz]
from out_mominxy import phiinxy
it = 250
iz = 48
phiinxy(it, iz, xr_phi, flag="display")


# In[ ]:


# ### Tips: To analyze interactively
# **1. xr_phi is just a xarray data set. `print(xr_phi)` shows the coordinate information.**

# In[2]:


print(xr_phi)


# **2. Please see the help of functions, which describes its feature, required inputs and outputs.**

# In[3]:


help(phiinxy)


# **3. `flag="display"` is convenient to see the result. On the other hand, choose `flag=None` to get data as a Numpy array and to draw figure as you like.**

# In[4]:


import numpy as np
import matplotlib.pyplot as plt
it = 250
iz = 48
data = phiinxy(it, iz, xr_phi, flag=None)
xx = data[:,:,0]
yy = data[:,:,1]
phixy = data[:,:,2]
plt.pcolormesh(xx,yy,phixy,shading="auto")
plt.xlim(-30,30)
plt.ylim(-30,30)
plt.show()


# **4. If you are familiar with Python, you can analyze everything by yourself. However, be careful that you also need to be familiar with the formulation of local
# flux-tube gyrokinetic simulation.**
# For example, flux-surface-averaged electric field spectrum is not $\int \frac{k_\perp^2}{2}|\phi_{k_x,ky}(z)|^2 dz/(\int dz)$, but $\int \sqrt{g} \frac{k_\perp^2}{2}|\phi_{k_x,ky}(z)|^2 dz/(\int \sqrt{g} dz)$ including the Jacobian of the employed flux-tube coordinates.
# The squared wavenumber is not $k_x^2 + k_y^2$, but $k_\perp^2 = g^{xx} k_x^2 + 2g^{xy}k_xk_y + g^{yy}k_y^2$, including the metrics.
#
# See GKV manual for theoretical discription, and use `diag_geom` module to pick up geometrical quantities.

# In[5]:


from diag_geom import ksq
from diag_intgrl import intgrl_thet

zz=xr_phi.zz
ky=xr_phi.ky
kx=xr_phi.kx
phi=xr_phi.phi[-1,:,:,:]
print(phi)

wrong_ksq = kx**2+ky**2
wrong_energy = np.sum(0.5*np.abs(phi)**2*wrong_ksq,axis=0)/len(zz)
proper_energy = intgrl_thet(0.5*ksq*np.abs(phi)**2)
plt.plot(ky,wrong_energy,"x-",label=r"Wrong energy")
plt.plot(ky,proper_energy,"o-",label=r"Proper energy")
plt.xlabel(r"$k_y$")
plt.ylabel(r"Energy spectrum $\int \sqrt{g} k_\perp^2|\phi_{k_x,ky}(z)|^2/2 dz/(\int \sqrt{g} dz)$")
plt.legend()
plt.show()


# In[ ]: