#!/usr/bin/env python # coding: utf-8 # # Ch05: Magnetohydrodynamics # # *Plasma as a fluid* # We now focus on the plasma as a single fluid. This mean averaging the electron motion over large spatial and time scales so the electron physics can be encapsulated readily into the ion flow. The trick is to retain as much as possible of the electron physics. But not too much, or the code will have to resolve electron effects which are much smaller and much faster than the ion fluid. # # This would reduce the numerical time step in our discrete equations, making our code very slow. This chapter unite both ion and electron physics via a generalized Ohm's law and the global electromagnetic fields, as done by Vlasov for kinetic theory. # The first step in getting the single fluid approximation is to sum up all the quantities evolved by the fluid equations to get # $$\begin{align} # &n=\sum_{i,e}n_{i,e}\tag{5.1}\\ # &\rho=\sum_{i,e}n_{i,e}m_{i,e}\tag{5.2}\\ # &n\vec u=\frac{1}{\rho}\sum_{i,e}\rho_{i,e}u_{i,e}\tag{5.3}\\ # &\vec j=Z_ien_i\vec u_i-en_e\vec u_e\tag{5.4}\end{align}$$ # From now on, we will call $\vec u$ the bulk flow of the fluid. # ## Generalized Ohm's law # We now go beyond the simple understanding that electric field can generate only electrical currents inside a plasma, or vise versa,i.e. # $$\vec E=\eta \vec j$$ # In fact, we have seen that there is only a relashionship between currents and electric field through resistivity. When the resistivity $\eta$ is simply 0, there is no electric field associated with electrical currents. But there are other cause for electric field to exists inside a plasma: # 1. Dynamo # 2. Hall effect # 3. Electron pressure (also known as the Biermann battery term) # 4. Electron inertia # # Let's look at each of them. # ### Dynamo and the Hall effect # This electric field is created by the drift of electron fluid compare to the ion fluid caused by a magnetic field. For one particle we know that $ \vec F=q\vec E+q\vec v\times\vec B$. if we suppose that the electron fluid is tied to the ion fluid then the global force on the electron fluid statistically balance the local force $\vec F$ so we get $q\vec E+q\vec u_e\times\vec B=\vec 0$ which give the first electric field to take into account inside plasma: # $$\vec E=-\vec u_e\times\vec B\tag{5.5}$$ # Now we work on the single fluid approximation of plasma so we need to get rid of $\vec u_e$. The global current density is given by $$\vec j=Z_ien_i\vec u_i-en_e\vec u_e\tag{5.6}$$ # Using quasi-neutrality, we can rewrite $\vec u_e$ as # $$u_e=\vec u_i-\frac{1}{en_e}\vec j$$ # Since the ion mass $m_i$ is much larger than the electron mass $m_e$ # $$\vec u_i=\frac{\rho\vec u+m_e\vec j/e}{\rho}\approx\vec u\tag{5.7}$$ # We end up with $$\vec E=-\vec u\times\vec B+\frac{1}{en_e}\vec j\times\vec B\tag{5.8}$$ # The first term on the right hand side the the dynamo term. A conductive flow moving inside a magnetic field creates an electric field, call the dynamo field, after Michael Faraday who coined the term for a machine capable of generating electric fields from flowing conductive fluid inside a magnetic field. The second term is the Hall effect, named after Edwin Hall, who discovered that electric field are generated when current flows inside a magnetic field. # Eq. $(5.8)$ shows that the electric field in the frame of reference of the flowing plasma is not $\vec E$. It is in fact $\vec E\,^*$, given by # $$\vec E\,^*=\vec E+\vec u\times\vec B\tag{5.9}$$ # So, this state of affair requires us to rewrite the Hall effect as # $$\vec E\,^*=\frac{1}{en_e}\vec j\times\vec B\tag{5.10}$$ # and, of course, Ohm's law as # $$\vec E\,^*=\eta\vec j\tag{5.11}$$ # ### Electron pressure (a.k.a. the Biermann battery) # Electric fields can also create a pressure gradient of electron. if there is more electron at one location inside the plasma than at another, then there is an electrostatic field building up. Normally, the electric field should be able to restore electrostatic balance inside the plasma. What happens if it can't? Typically this situation arises when the thermal energy of the electrons is large enough to fight the electric field. In this case, the electric force density $\vec f_E=-en_e\vec E\,^*$ is balanced exactly by the pressure gradient caused by the difference in electron density $\vec f_{\nabla p}=\vec \nabla p_e$. Since we are at equilibrium we get # $$\vec E\,^*=-\frac{1}{en_e}\vec \nabla p_e\tag{5.12}$$ # This situation can also happen when there is a spacial difference in electron temperature between two regions. The region where the temperature is larger pushes against the region where the temperature is smaller, again generating a pressure gradient. This gradient can exist only if there is an electric field that can balance it. # # # This effect is of particular interest in astrophysics since it can generate magnetic fields. Starting from the induction equation $$\frac{\partial \vec B}{\partial t}=-\vec \nabla\times\vec E\tag{5.13}$$ and using the fact that $p_e=en_eT_e$. Note that we used $e$ here because the temperature $T_e$ is in $V$. # we get $$\vec \nabla p_e=eT_e\vec \nabla n_e+en_e\vec\nabla T_e\tag{5.14}$$ # Using Eq. $(5.13)$ and $(5.14)$ we obtain the following equation # $$\frac{\partial \vec B}{\partial t}=\vec \nabla\times\Big(\frac{T_e}{n_e}\vec \nabla n_e+\vec\nabla T_e\Big)\tag{5.15}$$ # We know that $\vec\nabla\times(\vec\nabla \varphi)=\vec 0$ and $\vec\nabla\times(\varphi\vec A)=\varphi(\vec\nabla\times\vec A)+(\vec\nabla \varphi)\times\vec A$. So we can simplify Eq. $(5.15)$ # $$\frac{\partial \vec B}{\partial t}=\vec \nabla\Big(\frac{T_e}{n_e}\Big)\times\vec \nabla n_e$$ # # and rewrite it as # $$\frac{\partial \vec B}{\partial t}=\Big(\frac{\vec \nabla T_e}{n_e}-\frac{T_e\vec \nabla n_e}{n_e^2}\Big)\times\vec \nabla n_e$$ # # We finally get the the Biermann battery term # $$\frac{\partial \vec B}{\partial t}=\frac{1}{en_e}\vec \nabla (eT_e)\times\vec \nabla n_e\tag{5.16}$$ # So magnetic fields are generated from gradients in density and temperature with perpendicular components. # ### Electron Inertia # We mentioned earlier that the electric field seen by charged particles inside a flow is not $\vec E$ but $\vec E\,^*$. This electric field will generate a pull on all charged particles. But the electrons are much lighter than the ions. Consequently, the action on the electron is most noticeable. If we ignore effects on the ions, we get $$-e\vec E\,^*=m_e\frac{\partial(\vec v_e-\vec u)}{\partial t}.$$ We had to remove the speed $\vec u$ of the fluid here since we are using the electric field in the frame of reference of the plasma $\vec E\,^*$. We could rewrite this equation as # $$e\vec E\,^*=\frac{m_e}{en_e}\frac{\partial \vec j}{\partial t}\tag{5.17}.$$ # However, this approach is not completely correct. # To use $\vec j$ we have turned the speed of the individual electrons $\vec v_e$ into the flow speed of the electron flow $\vec u_e$. Yet, we would like Eq. $(5.17)$ to be an expression of conservation of momentum for the electron flow inside the bulk flow. So we can simply rewrite it as a conservation equation by adding the bulk flow convective term $\nabla(\vec j\vec u)$ (i.e. $\vec j$ is advected by $\vec u$). Yet, the flow of electron does also advect its own momentum and this is where it becomes slightly more complicated and we just give the final answer here # $$\vec E\,^*=\frac{m_e}{e^2n_e}\Big(\frac{\partial \vec j}{\partial t}+\vec\nabla\cdot\big[\vec u\,\vec j+\vec j\,\vec u+\frac{1}{en_e}\vec j\,\vec j\big]\Big)\tag{5.18}$$ # ## The equations of extended magnetohydrodynamics # Extended magnetohydrodynamics is a particular version of magnetohydrodynamics where only low frequency electron physics is kept. While it does not capture microscopic electron physics like plasma oscillations, it capture most of the electron physics that is expressed on the ion spatial and time scales. # # We can now gather all the electric field terms from Eqs. $(5.10)$, $(5.11)$, $(5.12)$ and $(5.18)$ together to get the total electric field $\vec E\,^*$ inside the bulk flow is # $$\vec E\,^*=\frac{1}{en_e}\vec j\times\vec B+\eta\vec j-\frac{1}{en_e}\vec \nabla p_e+\frac{m_e}{e^2n_e}\Big(\frac{\partial \vec j}{\partial t}+\vec\nabla\cdot\big[\vec u\,\vec j+\vec j\,\vec u+\frac{1}{en_e}\vec j\,\vec j\big]\Big)$$ # As a result, the electric field as seen by the plasma in the rest frame is # $$\vec E=-\vec u\times\vec B+\eta\vec j+\frac{1}{en_e}\big(\vec j\times\vec B-\vec \nabla p_e\big)+\frac{m_e}{e^2n_e}\Big(\frac{\partial \vec j}{\partial t}+\vec\nabla\cdot\big[\vec u\,\vec j+\vec j\,\vec u+\frac{1}{en_e}\vec j\,\vec j\big]\Big)\tag{5.19}$$ # This term corresponds to the time it takes for the electrons to respond to changes in electric fields due to their mass. # If you remember Ampere's laws # $$\frac{\partial\vec {E}}{\partial t}=\frac{1}{\mu_0\varepsilon_0}\vec \nabla\times\vec {B}-\frac{1}{\varepsilon_0}\vec {j}\tag{5.20}$$ # and Faraday's law # $$\frac{\partial\vec {B}}{\partial t}=-\vec \nabla\times\vec {E}\tag{5.21}$$ # from the previous chapter we might be able to see a little problem arising. # While it can be shown that $\vec E$ can be a solution of Eq. $(5.20)$ while satisfying the constraint of Maxwell's equations given in $(4.7)$, there is no mathematical proof that $\vec E$ can satisfy simultaneously Eq. $(5.19)$ and $(5.20)$. Well this might be problematic. Now we have two options. The first is to ignore radiations. So we can rewrite Eq. $(5.20)$ as # $$\vec \nabla\times\vec {B}=\mu_0\vec {j}\tag{5.22}$$ # In this case all is well as long as we enforce that $\vec\nabla\cdot \vec B=0$. # The other possibility is to write Eq. $(5.19)$ in its original form, that is # $$\frac{\partial \vec j}{\partial t}+\vec\nabla\cdot\big[\vec u\,\vec j+\vec j\,\vec u+\frac{1}{en_e}\vec j\,\vec j\big]=\frac{e^2n_e}{m_e}\Big(\vec E+\vec u\times\vec B-\eta\vec j-\frac{1}{en_e}[\,\vec j\times\vec B-\vec \nabla p_e]\Big)\tag{5.23}$$ # Now we have obtained a time advanced equation for the current density $\vec j$. Yet this equation is still on the electron time scale and, if solved directly, will suffer from a very small CFL conditions. However, it is possible to solve it implicitly when the value $e^2n_e/m_e$ is relatively large ($n_e>10^{10}m^{-3}$). In this case, we can force the right hand side of Eq. $(5.23)$ to be 0. # At this point, we can first solve Maxwell's equations using # $$\begin{align} # &\frac{\partial\vec {E}}{\partial t}=\frac{1}{\mu_0\varepsilon_0}\vec \nabla\times\vec {B}-\frac{1}{\varepsilon_0}\vec {j}\\ # &\frac{\partial\vec {B}}{\partial t}=-\vec \nabla\times\vec {E}\end{align}$$ # together with # $$\frac{\partial \vec j}{\partial t}+\vec\nabla\cdot\big[\vec u\,\vec j+\vec j\,\vec u+\frac{1}{en_e}\vec j\,\vec j\big]=0\tag{5.24}$$ # on the ion time scale. # Then, we correct the values of $\vec E$ and $\vec j$ using # $$\vec E+\vec u\times\vec B-\eta\vec j-\frac{1}{en_e}[\,\vec j\times\vec B-\vec \nabla p_e]=0\tag{5.25}$$ # At this point, we have implicitly found the contribution of electron physics *at the ion time scale* using the generalized Ohm's law, without the electron inertia term. One important property of Eq. $(5.25)$ is its lack of derivatives. The coupling between $E$ and $j$ forms a linear set of equations, which matrix is only $3\times3$ and can be inverted analytically. # Using the form of energy density defined in Ch04 # $$m_i\epsilon=\frac{3}{2}p+\frac{1}{2}\rho\vec u^2\tag{5.26}$$ # where $m_i$ is the ion mass, # our final set of equations for extended magnetohydrodynamics are # $$\begin{align} # &\frac{\partial n}{\partial t}+\vec \nabla . \big(n\vec u\big)=0\\ # &\frac{\partial n\vec u}{\partial t}+\vec \nabla . \big(n\vec u\vec u\big)=\vec j\times\vec B-\vec\nabla\cdot p\\ # &\frac{\partial \epsilon}{\partial t}+\vec \nabla \cdot (\epsilon\vec u)=\frac{1}{m}\big(-\vec \nabla \cdot Q+\vec E\cdot\vec j\big)\\ # &\frac{\partial\vec {E}}{\partial t}=\frac{1}{\mu_0\varepsilon_0}\vec \nabla\times\vec {B}-\frac{1}{\varepsilon_0}\vec {j}\\ # &\frac{\partial\vec {B}}{\partial t}=-\vec \nabla\times\vec {E}\\ # &\frac{\partial \vec j}{\partial t}+\vec\nabla\cdot\big[\vec u\,\vec j+\vec j\,\vec u+\frac{1}{en_e}\vec j\,\vec j\big]=0\\ # &\vec E+\vec u\times\vec B-\eta\vec j-\frac{1}{en_e}[\,\vec j\times\vec B-\vec \nabla p_e]=0 # \end{align}$$ # ## Resistive magnetohydrodynamics # If the value $e^2n_e/m_e$ is too small ($n_e<<10^{10}m^{-3}$) then the model presented above is not strictly correct. While it can still be used, it may not be possible to obtain a physical solution. If electron physics is important then, we need to use the two-fluid magnetohydrodynamics equations discussed in the previous chapter. However, if electron physics cna be neglected, then we can drop Eq. $(5.25)$ and simplify Eq. $(5.24)$ to a trimmed down version of the generalized Ohm's law # $$\vec E=-\vec u\times\vec B+\eta\vec j\tag{5.27}$$ # # Since we cannot that the electric field verifies both Maxwell's equations and Ohm's law, because we lost the degrees of freedom provided bu the time evolution of the current density, we need to make so further concessions and drop electromagnetic waves, i.e. $\partial\vec E/\partial t=\vec 0$. At this point the current density is fully defined by # $$\vec j=\frac{1}{\mu_0}\vec\nabla\times\vec B\tag{5.28}$$ # The other equations remain unchanged after we have replaced both $\vec E$ and $\vec j$ by their values given by Ohm's law and ampere's law. # $$\begin{align} # &\frac{\partial n}{\partial t}+\vec \nabla . \big(n\vec u\big)=0\\ # &\frac{\partial n\vec u}{\partial t}+\vec \nabla . \big(n\vec u\vec u\big)=\frac{1}{\mu_0}\big(\vec\nabla\times\vec B\big)\times\vec B-\vec\nabla\cdot p\\ # &\frac{\partial \epsilon}{\partial t}+\vec \nabla \cdot (\epsilon\vec u)=\frac{1}{m}\Big(-\vec \nabla \cdot Q+\frac{1}{\mu_0}\big(\vec\nabla\times\vec B\big)\cdot\big[\frac{\eta}{\mu_0}\vec\nabla\times\vec B-\vec u\times\vec B\big]\Big)\\ # &\frac{\partial\vec {B}}{\partial t}=-\vec \nabla\times\big(\frac{\eta}{\mu_0}\vec\nabla\times\vec B-\vec u\times\vec B\big)\\ # &\vec \nabla\cdot\vec B=0\\ # \end{align}$$ # Here both $\vec E$ and $\vec j$ have disappeared from the equations. # ## Ideal magnetohydrodynamics # Ideal magnetohydrodynamics suppose the resistivity $\eta$ negligible compared to dynamo. And we are left with # $$\begin{align} # &\frac{\partial n}{\partial t}+\vec \nabla . \big(n\vec u\big)=0\\ # &\frac{\partial n\vec u}{\partial t}+\vec \nabla . \big(n\vec u\vec u\big)=\frac{1}{\mu_0}\big(\vec\nabla\times\vec B\big)\times\vec B-\vec\nabla\cdot p\\ # &\frac{\partial \epsilon}{\partial t}+\vec \nabla \cdot (\epsilon\vec u)=-\frac{\vec \nabla \cdot Q}{m}\\ # &\frac{\partial\vec {B}}{\partial t}=\vec \nabla\times\big(\vec u\times\vec B\big)\\ # &\vec\nabla\cdot\vec B=0 # \end{align}$$ # ## Extended magnetohydrodynamics code # It is now possible to construct an extended magnetohydrodynamics code using the numerical method presented in Ch04 and the equations of extended magnetohydrodynamics. # ### Python housekeeping # In[1]: import numpy as np import math from math import sqrt import matplotlib import matplotlib.pyplot as plt import matplotlib.colors as colors from numba import jit from ipywidgets.widgets.interaction import show_inline_matplotlib_plots from IPython.display import clear_output from ipywidgets import FloatProgress from IPython.display import display # ### Definitions and constants # as before, we define some constants that will be use to remember which component in the array correspond to which quantity, e.g. the temperature can be accessed using `tp` at the location `i,j` as `Q[tp,i,j]`. `rh` is the number density, `mx` is the momentum density in the x-direction, `my` and `mz` in the y and z directions. `en` is for the energy density. We also introduce new variables: # 1. the magnetic field constants: `bx`, `by`, `bz` # 2. the electric field constants: `ex`, `ey` and `ez` # 3. the current density constants: `jx`, `jy` and `jz` # 4. the energy density for the electrons: `et` # 6. the electron density: `ne` # 7. the electron pressure: `ep` # 8. the extended magneto hydrodynamics resistivity: `etax` # 9. the electron and ion temperatures, `tpe` and `tpi` respectively # `nQ` is the total number of variables, from `rh` to `tp`. # In[2]: rh=0 mx=1 my=2 mz=3 en=4 bx=5 by=6 bz=7 ex=8 ey=9 ez=10 jx=11 jy=12 jz=13 et=14 ne=15 ep=16 etax=17 tpi=18 tpe=19 nQ=20 # `mu` is the atomic number of the fluid, in this case aluminum. `aindex` is the adiabatic index. `Z` is the average ionization factor and `er2` is used to reduce the speed of light, so computations run faster. While this is physically incorrect to reduce the speed of light, it does nto cause problems as long as the fastest speeds inside the plasmas are much smaller than the speed of light. # In[3]: Z=4. mu=27. er2=100. aindex=1.1 # Here we define our dimensional parameters. `L0` is the geometrical scale length, `t0` the time scale, `n0` the density scale, `v0` is the velocity scale, `p0` the pressure scale, `te0` the temperature scale. We also define a couple of floor quantities under which no density, temperature and pressure can go. Then we define all our electromagnetic related constants. # In[4]: L0=1.0e-3 t0=100.0e-9 n0=6.0e28 v0=L0/t0 p0=mu*1.67e-27*n0*v0**2 te0=p0/n0/1.6e-19 rh_floor = 1.0E-8 T_floor = 0.026/te0 rh_mult = 1.1 P_floor = rh_floor*T_floor mu0=4.e-7*np.pi eps0=8.85e-12 mime=mu*1836. memi=1./mime B0=math.sqrt(mu*1.67e-27*mu0*n0)*v0 E0=v0*B0 J0=B0/(mu0*L0) eta0=mu0*v0*L0 cflm = 0.8 gma2 = 2.8e-8*mu0*L0**2*n0 lil0 = 2.6e5*math.sqrt(mu/n0/mu0)/L0 cab = t0/L0/math.sqrt(eps0*mu0) ecyc = 1.76e11*B0*t0 clt = cab/er2 cfva2=(cab/er2)**2 vhcf = 1.0e6 sigma = 1.2e2 taui0 = 6.2e14*sqrt(mu)/n0 omi0 = 9.6e7*Z*B0/mu bapp = 0. visc0 = 5.3*t0 ome0 = mime*omi0 taue0 = 61.*taui0 # We now define our geometrical functions. They are used quite often so we will compile them inside the code to accelerate computations. To avoid writing two codes, one for slab symmetry and one for cylindrical symmetry, we wrote this code in for slab symmetry and in Cartesian coordinates, i.e. $(x,y,z)$. When using the code for problem with cylindrical symmetry then $$\begin{align} # &x\rightarrow r\\ # &y\rightarrow z\\ # &z\rightarrow \phi\end{align}$$ # # The `rc` function computes the radius in cylindrical coordinates, which appears in source terms and in the conservation equations. Clearly, this is simply $x$ as mentioned above. That's why `rc` returns `xc`. To turn the code into slab symmetry, just have the function `rc` return `1`. # In[5]: @jit(nopython=True) def xc(i): return loc_lxd + (i-ngu+1)/dxi @jit(nopython=True) def rc(i): return xc(i) @jit(nopython=True) def yc(j): return loc_lyd + (j-ngu+1)/dyi # We now add radius and angle function for axi-symmetric distribution in slab geometries, not ot be confused with `rc` above. Note the the letter `c`, for center, is now gone. # In[6]: @jit(nopython=True) def r(i,j): return math.sqrt((xc(i)-(lxd+lxu)/2)**2 + (yc(j)-(lyd+lyu)/2)**2) @jit(nopython=True) def rz(i,j): return sqrt(xc(i)**2+yc(j)**2) @jit(nopython=True) def theta(i,j): return math.atan(yc(j)-(lyd+lyu)/2,xc(i)-(lxd+lxu)/2) @jit(nopython=True) def Icur(t): return Ipeak*math.sin(np.pi*t/(2*tr))*math.exp(-t/tdp) # ### Finite volume method and time advance algorithms # We now compute the sources for all our conservation equations. # In[7]: def get_sources(Qin): sourcesin=np.zeros((nx,ny,nQ)) nuei=np.empty((nx,ny)) gyro = memi*ecyc linv = 1./lil0 fac = lil0/(mime + Z) Zi = 1./Z Zin = 1./(Z + 1.) dti = 0.2/dt dx2 = 1./(dt*dxi**2) dx = 1./dxi mui = (1. + Z)*t0 t0i = 1./t0 kap_e = 0. kap_i = 0. for j in range(ngu, ny-ngu): for i in range(ngu, nx-ngu): rbp = 0.5*(rc(i+1) + rc(i)) rbm = 0.5*(rc(i) + rc(i-1)) dni = Qin[i,j,rh] dne = Qin[i,j,ne] dnii = 1./dni dnei = 1./dne vix = Qin[i,j,mx]*dnii viy = Qin[i,j,my]*dnii viz = Qin[i,j,mz]*dnii vex = vix - lil0*Qin[i,j,jx]*dnei vey = viy - lil0*Qin[i,j,jy]*dnei vez = viz - lil0*Qin[i,j,jz]*dnei Qin[i,j,tpi] = (aindex - 1)*(Qin[i,j,en]*dnii - 0.5*(vix**2 + viy**2 + viz**2)) if(Qin[i,j,tpi] < T_floor): Qin[i,j,tpi] = T_floor Qin[i,j,tpe] = (aindex - 1)*(Qin[i,j,et]*dnei - 0.5*memi*(vex**2 + vey**2 + vez**2)) if(Qin[i,j,tpe] < T_floor): Qin[i,j,tpe] = T_floor eta_s = 0.5*Z/(Zin* Qin[i,j,tpi])**1.5 Qin[i,j,etax] = 1./sqrt(1./eta_s**2 + 1.e6*dni**2) nuei[i,j] = gma2*dne*Qin[i,j,etax] sourcesin[i,j,ne] = dti*(Z*Qin[i,j,rh] - Qin[i,j,ne]) sourcesin[i,j,mx] = Qin[i,j,jy]*Qin[i,j,bz] - Qin[i,j,jz]*Qin[i,j,by] sourcesin[i,j,my] = Qin[i,j,jz]*Qin[i,j,bx] - Qin[i,j,jx]*Qin[i,j,bz] sourcesin[i,j,mz] = Qin[i,j,jx]*Qin[i,j,by] - Qin[i,j,jy]*Qin[i,j,bx] ux = vix uy = viy uz = viz P = dni*Qin[i,j,tpi] + dne*Qin[i,j,tpe] sourcesin[i,j,en] = Z*gyro*dni*(vix*Qin[i,j,ex] + viy*Qin[i,j,ey] + viz*Qin[i,j,ez])- memi*lil0*nuei[i,j]*(vix*Qin[i,j,jx] + viy*Qin[i,j,jy] + viz*Qin[i,j,jz])+ 3*memi*nuei[i,j]*dni*(Qin[i,j,tpe] - Qin[i,j,tpi]) sourcesin[i,j,et] = -gyro*dne*(vex*Qin[i,j,ex] + vey*Qin[i,j,ey] + vez*Qin[i,j,ez])+ memi*lil0*nuei[i,j]*(vix*Qin[i,j,jx] + viy*Qin[i,j,jy] + viz*Qin[i,j,jz]) - 3*memi*nuei[i,j]*dni*(Qin[i,j,tpe] - Qin[i,j,tpi]) Qin[i,j,ep] = dne*Qin[i,j,tpe] sourcesin[i,j,ne] = sourcesin[i,j,ne]*rc(i) sourcesin[i,j,mx] = sourcesin[i,j,mx]*rc(i) + (Qin[i,j,mz]*viz + P) sourcesin[i,j,mz] = sourcesin[i,j,mz]*rc(i) - Qin[i,j,mz]*vix sourcesin[i,j,en] = sourcesin[i,j,en]*rc(i) sourcesin[i,j,et] = sourcesin[i,j,et]*rc(i) sourcesin[i,j,jx] = (uz*Qin[i,j,jz] + Qin[i,j,jz]*uz - lil0*dnei*Qin[i,j,jz]*Qin[i,j,jz]) sourcesin[i,j,jz] = -(ux*Qin[i,j,jz] + Qin[i,j,jx]*uz - lil0*dnei*Qin[i,j,jz]*Qin[i,j,jx]) sourcesin[i,j,bz] = Qin[i,j,ey] sourcesin[i,j,ez] = -cfva2*Qin[i,j,by] return sourcesin # We compute the electric field using constrained transport, to make sure that $\vec\nabla\cdot\vec B=0$ # In[8]: @jit(nopython=True) def get_ect(flux_x,flux_y): ect=np.zeros((nx,ny)) for j in range(ngu-1, ny-ngu+1): for i in range(ngu-1, nx-ngu): ect[i,j] = 0.25*(-flux_x[i-1,j,by] - flux_x[i,j,by] + flux_y[i,j-1,bx] + flux_y[i,j,bx]) #ect[:,1] = ect[:,2] #ect[1,:] = -ect[2,:] return ect # Here do the time advance of the conservation equations. # In[9]: @jit(nopython=True) def advance_time_level(Qin,ect,flux_x,flux_y,source): Qout=np.zeros((nx,ny,nQ)) fac = 0*mime/lil0 dxt = dxi*dt dyt = dyi*dt ect=get_ect(flux_x,flux_y) for j in range(ngu, ny-ngu): for i in range(ngu, nx-ngu): rbp = 0.5*(rc(i+1) + rc(i)) rbm = 0.5*(rc(i) + rc(i-1)) rci = 1./rc(i) Qout[i,j,rh] = Qin[i,j,rh]*rc(i) - dxt*(flux_x[i,j,rh]*rbp-flux_x[i-1,j,rh]*rbm) - dyt*(flux_y[i,j,rh]-flux_y[i,j-1,rh])*rc(i) Qout[i,j,mx] = Qin[i,j,mx]*rc(i) - dxt*(flux_x[i,j,mx]*rbp-flux_x[i-1,j,mx]*rbm) - dyt*(flux_y[i,j,mx]-flux_y[i,j-1,mx])*rc(i) + dt*source[i,j,mx] Qout[i,j,my] = Qin[i,j,my]*rc(i) - dxt*(flux_x[i,j,my]*rbp-flux_x[i-1,j,my]*rbm) - dyt*(flux_y[i,j,my]-flux_y[i,j-1,my])*rc(i) + dt*source[i,j,my] Qout[i,j,mz] = Qin[i,j,mz]*rc(i) - dxt*(flux_x[i,j,mz]*rbp-flux_x[i-1,j,mz]*rbm) - dyt*(flux_y[i,j,mz]-flux_y[i,j-1,mz])*rc(i) + dt*source[i,j,mz] Qout[i,j,en] = Qin[i,j,en]*rc(i) - dxt*(flux_x[i,j,en]*rbp-flux_x[i-1,j,en]*rbm) - dyt*(flux_y[i,j,en]-flux_y[i,j-1,en])*rc(i) + dt*source[i,j,en] Qout[i,j,bx] = Qin[i,j,bx]*rc(i) - 0.5*dyt*(ect[i,j+1]*rc(i) - ect[i,j-1]*rc(i)) Qout[i,j,by] = Qin[i,j,by]*rc(i) + 0.5*dxt*(ect[i+1,j]*rc(i+1) - ect[i-1,j]*rc(i-1)) Qout[i,j,bz] = Qin[i,j,bz]*rc(i) - dxt*(flux_x[i,j,bz]*rbp-flux_x[i-1,j,bz]*rbm) - dyt*(flux_y[i,j,bz]-flux_y[i,j-1,bz])*rc(i) + dt*source[i,j,bz] Qout[i,j,ex] = Qin[i,j,ex]*rc(i) - dxt*(flux_x[i,j,ex]*rbp-flux_x[i-1,j,ex]*rbm) - dyt*(flux_y[i,j,ex]-flux_y[i,j-1,ex])*rc(i) + dt*source[i,j,ex] Qout[i,j,ey] = Qin[i,j,ey]*rc(i) - dxt*(flux_x[i,j,ey]*rbp-flux_x[i-1,j,ey]*rbm) - dyt*(flux_y[i,j,ey]-flux_y[i,j-1,ey])*rc(i) + dt*source[i,j,ey] Qout[i,j,ez] = Qin[i,j,ez]*rc(i) - dxt*(flux_x[i,j,ez]*rbp-flux_x[i-1,j,ez]*rbm) - dyt*(flux_y[i,j,ez]-flux_y[i,j-1,ez])*rc(i) + dt*source[i,j,ez] Qout[i,j,ne] = Qin[i,j,ne]*rc(i) - dxt*(flux_x[i,j,ne]*rbp-flux_x[i-1,j,ne]*rbm) - dyt*(flux_y[i,j,ne]-flux_y[i,j-1,ne])*rc(i) + dt*source[i,j,ne] Qout[i,j,jx] = Qin[i,j,jx]*rc(i) - dxt*(flux_x[i,j,jx]*rbp-flux_x[i-1,j,jx]*rbm) - dyt*(flux_y[i,j,jx]-flux_y[i,j-1,jx])*rc(i) + dt*source[i,j,jx] + fac*dxt*(flux_x[i,j,ep]-flux_x[i-1,j,ep])*rc(i) Qout[i,j,jy] = Qin[i,j,jy]*rc(i) - dxt*(flux_x[i,j,jy]*rbp-flux_x[i-1,j,jy]*rbm) - dyt*(flux_y[i,j,jy]-flux_y[i,j-1,jy])*rc(i) + dt*source[i,j,jy] + fac*dyt*(flux_y[i,j,ep]-flux_y[i,j-1,ep])*rc(i) Qout[i,j,jz] = Qin[i,j,jz]*rc(i) - dxt*(flux_x[i,j,jz]*rbp-flux_x[i-1,j,jz]*rbm) - dyt*(flux_y[i,j,jz]-flux_y[i,j-1,jz])*rc(i) + dt*source[i,j,jz] Qout[i,j,et] = Qin[i,j,et]*rc(i) - dxt*(flux_x[i,j,et]*rbp-flux_x[i-1,j,et]*rbm) - dyt*(flux_y[i,j,et]-flux_y[i,j-1,et])*rc(i) + dt*source[i,j,et] Qout[i,j,rh:ne+1] = Qout[i,j,rh:ne+1]*rci return Qout # After advancing the current using Eq. $(5.24)$ we compute $\vec E$ and $\vec j$ from Eq. $(5.25)$ # In[10]: @jit(nopython=True) def implicit_source2(Qin,flux_x,flux_y): fac = 1./(mime + Z) Ainv=np.zeros((3,3)) mb = ecyc*dt mb2 = mb**2 for j in range(ngu, ny-ngu): for i in range(ngu, nx-ngu): dne = Qin[i,j,ne] zero = 1. if(dne < Z*rh_floor*rh_mult): zero = 0. dnii = 1./Qin[i,j,rh] ux = Qin[i,j,mx]*dnii uy = Qin[i,j,my]*dnii uz = Qin[i,j,mz]*dnii hx = Qin[i,j,bx] hy = Qin[i,j,by] hz = Qin[i,j,bz] ma = 1. + dt*gma2*dne*(Qin[i,j,etax] + dt*cfva2) ma2 = ma**2 mamb = ma*mb denom = zero/(ma*(ma2 + mb2*(hx**2 + hy**2 + hz**2))) Ainv[0,0] = denom*(ma2 + mb2*hx**2) Ainv[1,1] = denom*(ma2 + mb2*hy**2) Ainv[2,2] = denom*(ma2 + mb2*hz**2) Ainv[0,1] = denom*(mb2*hx*hy - mamb*hz) Ainv[1,2] = denom*(mb2*hy*hz - mamb*hx) Ainv[2,0] = denom*(mb2*hz*hx - mamb*hy) Ainv[1,0] = denom*(mb2*hx*hy + mamb*hz) Ainv[2,1] = denom*(mb2*hy*hz + mamb*hx) Ainv[0,2] = denom*(mb2*hz*hx + mamb*hy) Qjx = Qin[i,j,jx] + dt*gma2*(dne*(Qin[i,j,ex] + uy*hz - uz*hy) + dxi*lil0*(flux_x[i,j,ep] - flux_x[i-1,j,ep])) Qjy = Qin[i,j,jy] + dt*gma2*(dne*(Qin[i,j,ey] + uz*hx - ux*hz) + dyi*lil0*(flux_y[i,j,ep] - flux_y[i,j-1,ep])) Qjz = Qin[i,j,jz] + dt*gma2*(dne*(Qin[i,j,ez] + ux*hy - uy*hx)) Qin[i,j,jx] = Ainv[0,0]*Qjx + Ainv[0,1]*Qjy + Ainv[0,2]*Qjz Qin[i,j,jy] = Ainv[1,0]*Qjx + Ainv[1,1]*Qjy + Ainv[1,2]*Qjz Qin[i,j,jz] = Ainv[2,0]*Qjx + Ainv[2,1]*Qjy + Ainv[2,2]*Qjz Qin[i,j,ex:ez] -= dt*cfva2*Qin[i,j,jx:jz] # We now write the function that computes the flux `flx` and the freezing speeds `cfx` along the x-direction. These speeds correspond to the largest speed found in the hyperbolic equation at the given location $(x,y)$. # In[11]: @jit(nopython=True) def calc_flux_x(Qx): cfx=np.zeros((nx,nQ)) flx=np.zeros((nx,nQ)) fac = 1./(mime + Z) for i in range(0, nx): dni = Qx[i,rh] #dnii =0 #if (dni>0): dnii = 1./dni vx = Qx[i,mx]*dnii vy = Qx[i,my]*dnii vz = Qx[i,mz]*dnii dne = Qx[i,ne] #dnei =0 #if (dne>0): dnei = 1./dne vex = vx - lil0*Qx[i,jx]*dnei vey = vy - lil0*Qx[i,jy]*dnei vez = vz - lil0*Qx[i,jz]*dnei hx = Qx[i,bx] hy = Qx[i,by] hz = Qx[i,bz] Pri = (aindex - 1)*(Qx[i,en] - 0.5*dni*(vx**2 + vy**2 + vz**2)) if(Pri < P_floor): Pri = P_floor Pre = (aindex - 1)*(Qx[i,et] - 0.5*memi*dne*(vex**2 + vey**2 + vez**2)) if(Pre < Z*P_floor): Pre = Z*P_floor P = Pri + Pre vx = (mime*vx + Z*vex)*fac vy = (mime*vy + Z*vey)*fac vz = (mime*vz + Z*vez)*fac flx[i,rh] = Qx[i,mx] flx[i,mx] = Qx[i,mx]*vx + P flx[i,my] = Qx[i,my]*vx flx[i,mz] = Qx[i,mz]*vx flx[i,en] = (Qx[i,en] + Pri)*vx flx[i,bx] = 0.0 flx[i,by] = -Qx[i,ez] flx[i,bz] = Qx[i,ey] flx[i,ex] = 0. flx[i,ey] = cfva2*hz flx[i,ez] = -cfva2*hy flx[i,ne] = dne*vex flx[i,jx] = vx*Qx[i,jx] + Qx[i,jx]*vx - lil0*dnei*Qx[i,jx]*Qx[i,jx] flx[i,jy] = vx*Qx[i,jy] + Qx[i,jx]*vy - lil0*dnei*Qx[i,jx]*Qx[i,jy] flx[i,jz] = vx*Qx[i,jz] + Qx[i,jx]*vz - lil0*dnei*Qx[i,jx]*Qx[i,jz] flx[i,et] = (Qx[i,et] + Pre)*vex flx[i,ep] = Pre asqr = sqrt(aindex*P*dnii) # vf1 = abs(vx) + asqr vf1 = sqrt(vx**2 + hz**2*dnii + aindex*P*dnii) vf2 = clt vf3 = abs(vex) + asqr cfx[i,rh] = vf1 cfx[i,mx] = vf1 cfx[i,my] = vf1 cfx[i,mz] = vf1 cfx[i,en] = vf1 cfx[i,ne] = vf3 cfx[i,jx] = vf3 cfx[i,jy] = vf3 cfx[i,jz] = vf3 cfx[i,et] = vf3 cfx[i,bx] = vf2 cfx[i,by] = vf2 cfx[i,bz] = vf2 cfx[i,ex] = vf2 cfx[i,ey] = vf2 cfx[i,ez] = vf2 return cfx,flx # We now write the function that computes the flux `fly` and the freezing speeds `cfy` along the y-direction. # In[12]: @jit(nopython=True) def calc_flux_y(Qy): cfy=np.zeros((ny,nQ)) fly=np.zeros((ny,nQ)) fac = 1./(mime + Z) for j in range(0, ny): dni = Qy[j,rh] #dnii =0 #if (dni>0): dnii = 1./dni vx = Qy[j,mx]*dnii vy = Qy[j,my]*dnii vz = Qy[j,mz]*dnii dne = Qy[j,ne] #dnei =0 #if (dne>0): dnei = 1./dne vex = vx - lil0*Qy[j,jx]*dnei vey = vy - lil0*Qy[j,jy]*dnei vez = vz - lil0*Qy[j,jz]*dnei hx = Qy[j,bx] hy = Qy[j,by] hz = Qy[j,bz] Pri = (aindex - 1)*(Qy[j,en] - 0.5*dni*(vx**2 + vy**2 + vz**2)) if(Pri < P_floor): Pri = P_floor Pre = (aindex - 1)*(Qy[j,et] - 0.5*memi*dne*(vex**2 + vey**2 + vez**2)) if(Pre < Z*P_floor): Pre = Z*P_floor P = Pri + Pre vx = (mime*vx + Z*vex)*fac vy = (mime*vy + Z*vey)*fac vz = (mime*vz + Z*vez)*fac ## end if fly[j,rh] = Qy[j,my] fly[j,mx] = Qy[j,mx]*vy fly[j,my] = Qy[j,my]*vy + P fly[j,mz] = Qy[j,mz]*vy fly[j,en] = (Qy[j,en] + Pri)*vy fly[j,bx] = Qy[j,ez] fly[j,by] = 0.0 fly[j,bz] = -Qy[j,ex] fly[j,ex] = -cfva2*hz fly[j,ey] = 0. fly[j,ez] = cfva2*hx fly[j,ne] = dne*vey fly[j,jx] = vy*Qy[j,jx] + Qy[j,jy]*vx - lil0*dnei*Qy[j,jy]*Qy[j,jx] fly[j,jy] = vy*Qy[j,jy] + Qy[j,jy]*vy - lil0*dnei*Qy[j,jy]*Qy[j,jy] fly[j,jz] = vy*Qy[j,jz] + Qy[j,jy]*vz - lil0*dnei*Qy[j,jy]*Qy[j,jz] fly[j,et] = (Qy[j,et] + Pre)*vey fly[j,ep] = Pre asqr = sqrt(aindex*P*dnii) # vf1 = abs(vy) + asqr vf1 = sqrt(vy**2 + hz**2*dnii + aindex*P*dnii) vf2 = clt vf3 = abs(vey) + asqr cfy[j,rh] = vf1 cfy[j,mx] = vf1 cfy[j,my] = vf1 cfy[j,mz] = vf1 cfy[j,en] = vf1 cfy[j,ne] = vf3 cfy[j,jx] = vf3 cfy[j,jy] = vf3 cfy[j,jz] = vf3 cfy[j,et] = vf3 cfy[j,bx] = vf2 cfy[j,by] = vf2 cfy[j,bz] = vf2 cfy[j,ex] = vf2 cfy[j,ey] = vf2 cfy[j,ez] = vf2 return cfy,fly # # end def calc_flux_y # We now smooth out the different fluxes using a total diminishing variation scheme as done in Ch04 but for all variables # In[13]: def tvd2(Qin,n,ff,cfr): sl = 1 #use 0.75 for first order wr = cfr*Qin + ff wl = cfr*Qin - ff fr = wr fl = np.roll(wl,-1,axis=0) dfrp = np.roll(fr,-1,axis=0) - fr dfrm = fr - np.roll(fr,+1,axis=0) dflp = fl - np.roll(fl,-1,axis=0) dflm = np.roll(fl,+1,axis=0) - fl dfr=np.zeros((n,nQ)) dfl=np.zeros((n,nQ)) flux2 = tvd2_cycle(nQ, n, dfrp, dfrm, dflp, dflm, dfr, dfl, fr, fl, sl) return flux2 @jit(nopython=True) def tvd2_cycle(nQ, n, dfrp, dfrm, dflp, dflm, dfr, dfl, fr, fl, sl): flux2=np.zeros((n,nQ)) for l in range(nQ): for i in range(n): if(dfrp[i,l]*dfrm[i,l] > 0) : dfr[i,l] = dfrp[i,l]*dfrm[i,l]/(dfrp[i,l] + dfrm[i,l]) else: dfr[i,l] = 0 if(dflp[i,l]*dflm[i,l] > 0) : dfl[i,l] = dflp[i,l]*dflm[i,l]/(dflp[i,l] + dflm[i,l]) else: dfl[i,l] = 0 flux2 = 0.5*(fr - fl + sl*(dfr - dfl)) return flux2 # We truly compute all the fluxes here. # In[14]: def get_flux(Qin): flux_x=np.zeros((nx,ny,nQ)) flux_y=np.zeros((nx,ny,nQ)) fx=np.zeros((nx,nQ)) cfx=np.zeros((nx,nQ)) ffx=np.zeros((nx,nQ)) fy=np.zeros((ny,nQ)) cfy=np.zeros((ny,nQ)) ffy=np.zeros((ny,nQ)) for j in range(0, ny): cfx,ffx=calc_flux_x(Qin[:,j,:]) flux_x[:,j,:]=tvd2(Qin[:,j,:],nx,ffx,cfx) for i in range(0, nx): cfy,ffy=calc_flux_y( Qin[i,:,:]) flux_y[i,:,:]=tvd2(Qin[i,:,:],ny,ffy,cfy) return flux_x,flux_y # We now limit flows to avoid minor instabilities caused by round errors. # In[15]: def limit_flow(Qin): en_floor = rh_floor*T_floor/(aindex-1) for j in range(ngu-1, ny-ngu+1): for i in range(ngu-1, nx-ngu+1): if(Qin[i,j,rh] <= rh_floor or Qin[i,j,ne] <= Z*rh_floor): Qin[i,j,rh] = rh_floor Qin[i,j,ne] = Z*rh_floor Qin[i,j,mx] = 0.0 Qin[i,j,my] = 0.0 Qin[i,j,mz] = 0.0 Qin[i,j,en] = en_floor Qin[i,j,et] = Z*en_floor if(abs(Qin[i,j,jx]) > vhcf*Qin[i,j,rh]): sign=1 if (Qin[i,j,jx]<0): sign=-1; Qin[i,j,jx] = vhcf*Qin[i,j,rh]*sign if(abs(Qin[i,j,jy]) > vhcf*Qin[i,j,rh]): sign=1 if (Qin[i,j,jy]<0): sign=-1; Qin[i,j,jy] = vhcf*Qin[i,j,rh]*sign if(abs(Qin[i,j,jz]) > vhcf*Qin[i,j,rh]): sign=1 if (Qin[i,j,jz]<0): sign=-1; Qin[i,j,jz] = vhcf*Qin[i,j,rh]*sign if(abs(Qin[i,j,mx]) > clt*Qin[i,j,rh]): sign=1 if (Qin[i,j,mx]<0): sign=-1; Qin[i,j,mx] = clt*Qin[i,j,rh]*sign if(abs(Qin[i,j,my]) > clt*Qin[i,j,rh]): sign=1 if (Qin[i,j,my]<0): sign=-1; Qin[i,j,my] = clt*Qin[i,j,rh]*sign if(abs(Qin[i,j,mz]) > clt*Qin[i,j,rh]): sign=1 if (Qin[i,j,mz]<0): sign=-1; Qin[i,j,mz] = clt*Qin[i,j,rh]*sign if(MMask[i,j]==1): Qin[i,j,mx] = 0.0 Qin[i,j,my] = 0.0 Qin[i,j,mz] = 0.0 Qin[i,j,en] = en_floor Qin[i,j,et] = Z*en_floor if(BMask[i,j]==1): Qin[i,j,bz] = 0.0 # ### Domain and variable definitions # Everything you find below should be the only lines of code you need to change. # We now initialize the computational domain. `ngu` correspond to the number of points at the boundary. Leave it alone. `nx` and `ny` are the number of grid points along the x- and y-directions. `lxu` and `lxd` are the upper and lower corners of the mesh along the x-direction. `lyu` and `lyd` are the corners for the y-direction. `dxi` and `dyi` are the inverse of the grid resolution `dx` and `dy`. Since we use both `1/dx` and `/dy` a lot and it is more expansive to compute inverse, we compute them once and that's it. Since python is interactive, we could see them as constants here. # # **Every time you want to restart a new computation you need to recompile what's below.** # In[16]: ngu=2 nx=35 ny=40 lxu=15e-3/L0 lxd=0 lyd=0 lyu = lyd + (ny-2*ngu)*(lxu-lxd)/(nx-2*ngu) dxi = (nx-2*ngu)/(lxu-lxd) dyi = (ny-2*ngu)/(lyu-lyd) loc_lxd = lxd loc_lyd = lyd box=np.array([lxd,lxu,lyd,lyu])*L0/1e-3 # We now compute the maximum size of the time step using the $CFL$ condition # In[17]: def get_min_dt(Q): valf = 0. vmag = 0. for j in range(ny): for i in range(nx): valf0 = sqrt((Q[i,j,bx]**2 + Q[i,j,by]**2 + Q[i,j,bz]**2)/Q[i,j,rh]) if(valf0 > valf and Q[i,j,rh] > rh_mult*rh_floor): valf = valf0 vmag0 = sqrt((Q[i,j,mx]**2 + Q[i,j,my]**2))/Q[i,j,rh] if(vmag0 > vmag and Q[i,j,rh] > rh_mult*rh_floor): vmag = vmag0 cfl = 0.5*(1. - 0.6*(valf/clt)**2/(1. + (valf/clt)**2)) vmax = max(clt,vmag) dt_min = cfl/(dxi*vmax) return dt_min # Below we defined the initial `ti` and final `tf` times. `n_out` corresponds to to the total number of outputs. The rest are variable initialization. Note the most important one, the time `t`. # In[18]: ti = 00.0E-2 tf = 150.e-9/t0 tr = 150e-9/t0 tdp = 6*tr t = ti dt = cflm/(2.828*dxi*clt) n_steps=0 n_out=10 t_count=0 Ipeak=1e6/(J0*L0*L0) # And now all our arrays are defined here. Note that we use a mask `MMask` for parts of the simulation space were we want to preclude motion. That will come handy later in the course. # In[19]: Q=np.zeros((nx,ny,nQ)) MMask=np.zeros((nx,ny)) BMask=np.zeros((nx,ny)) Q1=np.copy(Q) Q2=np.copy(Q) sources=np.zeros((nx,ny,nQ)) flux_x=np.zeros((nx,ny,nQ)) flux_y=np.zeros((nx,ny,nQ)) ect=np.zeros((nx,ny)) # ### Initial conditions # That's for an empty domain. # In[20]: Q[:,:,rh]=rh_floor Q[:,:,en]=0.026/te0*Q[:,:,rh]*(aindex-1) Q[:,:,ne]=Z*Q[:,:,rh] Q[:,:,et]=Z*Q[:,:,en] K_rad=2e-3/L0 A_rad=10e-3/L0 height =2e-3/L0 for i in range (nx): for j in range (ny): if (yc(j)A_rad): MMask[i,j]=1 BMask[i,j]=1 Q[i,j,rh]=6e28/n0 if (yc(j)>=height and yc(j+1)K_rad and xc(i)(tf-ti)/n_out)): #everything below is plotting.... fig=plt.figure(figsize=(20, 19)) for i in range(1, rows+1): for j in range(1, columns+1): k=j+(i-1)*columns fig.add_subplot(rows, columns, k) if (k==1): data= np.flip(np.transpose(Q[ngu:nx-ngu-1,ngu:ny-ngu-1,rh]*n0),0) im = plt.imshow(data, cmap=plt.cm.nipy_spectral, norm=colors.LogNorm(vmin=data.min(), vmax=data.max()),extent=box) cb = fig.colorbar(im,fraction=0.046, pad=0.04) cb.ax.set_ylabel('$\log_{10}n_i(m^{-3})$', rotation=-90) if (k==3): data=Q[ngu:nx-ngu-1,ngu:ny-ngu-1,mx]**2 data+=Q[ngu:nx-ngu-1,ngu:ny-ngu-1,my]**2 data+=Q[ngu:nx-ngu-1,ngu:ny-ngu-1,mz]**2 data/=Q[ngu:nx-ngu-1,ngu:ny-ngu-1,rh]**2 data=np.sqrt(data)*L0/t0 data= np.flip(np.transpose(data),0) im = plt.imshow(data, extent=box) cb = fig.colorbar(im,fraction=0.046, pad=0.04) cb.ax.set_ylabel('$u(m/s)$', rotation=-90) if (k==2): data=Q[ngu:nx-ngu-1,ngu:ny-ngu-1,mx]**2 data/=Q[ngu:nx-ngu-1,ngu:ny-ngu-1,rh]**2 data=np.sqrt(data)*L0/t0 data= np.flip(np.transpose(data),0) im = plt.imshow(data, extent=box) cb = fig.colorbar(im,fraction=0.046, pad=0.04) cb.ax.set_ylabel('$u_r(m/s)$', rotation=-90) if (k==4): data=Q[ngu:nx-ngu-1,ngu:ny-ngu-1,my]**2 data/=Q[ngu:nx-ngu-1,ngu:ny-ngu-1,rh]**2 data=np.sqrt(data)*L0/t0 data= np.flip(np.transpose(data),0) im = plt.imshow(data, extent=box) cb = fig.colorbar(im,fraction=0.046, pad=0.04) cb.ax.set_ylabel('$u_z(m/s)$', rotation=-90) if (k==6): data=Q[ngu:nx-ngu-1,ngu:ny-ngu-1,mz]**2 data/=Q[ngu:nx-ngu-1,ngu:ny-ngu-1,rh]**2 data=np.sqrt(data)*L0/t0 data= np.flip(np.transpose(data),0) im = plt.imshow(data, extent=box) cb = fig.colorbar(im,fraction=0.046, pad=0.04) cb.ax.set_ylabel('$u_\phi(m/s)$', rotation=-90) if (k==5): data=Q[ngu:nx-ngu-1,ngu:ny-ngu-1,tpi]*te0 data= np.flip(np.transpose(data),0) im = plt.imshow(data, cmap=plt.cm.jet, extent=box) cb = fig.colorbar(im,fraction=0.046, pad=0.04) cb.ax.set_ylabel('$T(eV)$', rotation=-90) if (k==7): data=Q[ngu:nx-ngu-1,ngu:ny-ngu-1,jx]**2+Q[ngu:nx-ngu-1,ngu:ny-ngu-1,jy]**2+Q[ngu:nx-ngu-1,ngu:ny-ngu-1,jz]**2 data=np.sqrt(data)*J0 data= np.flip(np.transpose(data),0) im = plt.imshow(data, extent=box) cb = fig.colorbar(im,fraction=0.046, pad=0.04) cb.ax.set_ylabel('$J(A/m^2)$', rotation=-90) if (k==8): data=Q[ngu:nx-ngu-1,ngu:ny-ngu-1,bx]**2+Q[ngu:nx-ngu-1,ngu:ny-ngu-1,by]**2+Q[ngu:nx-ngu-1,ngu:ny-ngu-1,bz]**2 data=np.sqrt(data)*B0 data= np.flip(np.transpose(data),0) im = plt.imshow(data, extent=box) cb = fig.colorbar(im,fraction=0.046, pad=0.04) cb.ax.set_ylabel('$B(T)$', rotation=-90) im.set_interpolation('quadric') cb.ax.yaxis.set_label_coords(9, 0.5) plt.gca().axes.get_yaxis().set_visible(False) plt.gca().axes.get_xaxis().set_visible(False) if (j==1): plt.ylabel('z(mm)', rotation=90) plt.gca().axes.get_yaxis().set_visible(True) if (i==rows): plt.xlabel('r(mm)', rotation=0) plt.gca().axes.get_xaxis().set_visible(True) plt.tight_layout() clear_output() output_bar.value=0 display(progress_bar) # display the bar display(output_bar) # display the bar t_count=0 show_inline_matplotlib_plots() t=t+dt #we are now done. Let's turn off the progress bars. output_bar.close() del output_bar progress_bar.close() del progress_bar print("DONE") # © Pierre Gourdain, The University of Rochester, Charlie Seyler, Cornell University, 2020 # In[ ]: