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.
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:
Let's look at each of them.
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}$$
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.
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}$$
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}$$
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 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}$$
It is now possible to construct an extended magnetohydrodynamics code using the numerical method presented in Ch04 and the equations of extended magnetohydrodynamics.
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
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:
bx, by, bzex, ey and ezjx, jy and jzetneepetaxtpe and tpi respectivelynQ is the total number of variables, from rh to tp.
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.
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.
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.
@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.
@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)
We now compute the sources for all our conservation equations.
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$
@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.
@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)$
@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)$.
@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.
@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
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.
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.
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
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.
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
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.
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.
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))
That's for an empty domain.
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)<height):
if (xc(i)<K_rad or xc(i)>A_rad):
MMask[i,j]=1
BMask[i,j]=1
Q[i,j,rh]=6e28/n0
if (yc(j)>=height and yc(j+1)<height+2./dxi):
Q[i,j,rh]=6e25/n0
Q[i,j,en]=0.026/te0*Q[i,j,rh]*(aindex-1)
Q[i,j,ne]=Z*Q[i,j,rh]
Q[i,j,et]=Z*Q[i,j,en]
Below we define our boundary conditions. In this example, we do not use the time t and focus only on steady state solutions. The first block define the default boundary conditions for the whole domain. Note that the conditions for the y-direction are pretty simple. These boundary conditions correspond to an open boundary. everything can go through it.
Boundary conditions are a bit more complex for the x-direction. We are using a cylindrical coordinate system with open boundary conditions on the right but reflecting boundary conditions on the left. That's because nothing can pass through the axis of symmetry $r=0$.
At the end we have the boundary conditions that make up our problem, e.g. material free falling toward the axis of our coordinate system. Note that we are using very large density because the author works in high energy density plasmas. the initial parameters L0, t0 and n0 can be changed to accommodate other regimes. But, chances are that the CFL scaler cflm may have to be adjusted.
def set_bc(Qin,t):
Qin[nx-ngu,:,:] = Qin[nx-ngu-1,:,:]*xc(nx-ngu-1)/xc(nx-ngu)
Qin[nx-ngu+1,:,:] = Qin[nx-ngu,:,:]*xc(nx-ngu)/xc(nx-ngu+1)
Qin[:,1,:] = Qin[:,2,:]
Qin[:,0,:] = Qin[:,1,:]
Qin[:,ny-2,:] = Qin[:,ny-3,:]
Qin[:,ny-1,:] = Qin[:,ny-2,:]
for l in range(ngu):
for k in range(nQ):
if(k == mx or k == jx or k == ex or k == bx):
Qin[l,:,k] = -Qin[2*ngu-l-1,:,k]
if (k == mz or k == jz or k == ez or k == bz) :
Qin[l,:,k] = -Qin[2*ngu-l-1,:,k]
if (k == my or k == jy or k == ey or k == by) :
Qin[l,:,k] = Qin[2*ngu-l-1,:,k]
if (k == rh or k == ne or k == en or k == et) :
Qin[l,:,k] = Qin[2*ngu-l-1,:,k]
#our boundary condtions for injecting material radially into our domain
#note that all the values used are dimensionless
N=3
for i in range (nx):
for j in range (ngu):
if (xc(i)>K_rad and xc(i)<A_rad):
Qin[i,j,bz]=Icur(t)/(2*np.pi*xc(nx-1-k))
Let's define some plotting parameters
matplotlib.rcParams.update({'font.size': 22})
xi = np.linspace(lxd*L0, lxu*L0, nx-2*ngu-1)
yi = np.linspace(lyd*L0, lyu*L0, ny-2*ngu-1)
progress_bar = FloatProgress(min=ti, max=tf,description='Progress')
output_bar = FloatProgress(min=0, max=(tf-ti)/n_out,description='Next output')
columns = 2
rows = 4
box=np.array([lxd,lxu,lyd,lyu])*L0/1e-3
Here we go again!
xi = np.linspace(lxd*L0, lxu*L0, nx-2*ngu-1)
yi = np.linspace(lyd*L0, lyu*L0, ny-2*ngu-1)
while t<tf:
dt=get_min_dt(Q)
set_bc(Q,t)
sources=get_sources(Q)
flux_x,flux_y=get_flux(Q)
Q1=advance_time_level(Q,ect,flux_x,flux_y,sources)
implicit_source2(Q1,flux_x,flux_y)
limit_flow(Q1)
set_bc(Q1,t+dt)
sources=get_sources(Q1)
flux_x,flux_y=get_flux(Q1)
Q2=advance_time_level(Q1,ect,flux_x,flux_y,sources)
implicit_source2(Q2,flux_x,flux_y)
limit_flow(Q2)
Q=0.5*(Q+Q2)
limit_flow(Q)
n_steps+=1 #time step increased by 1
t_count+=dt #and the time by dt
progress_bar.value=t #we have defined some progress bars cause python is slow
output_bar.value+=dt # we need to monitor our progress and reduce resolution if it takes too long
if ((t==ti) or (t_count>(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")
FloatProgress(value=0.3011708828279219, description='Progress', max=1.5)
FloatProgress(value=0.0, description='Next output', max=0.15)