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
  5. the electron density: ne
  6. the electron pressure: ep
  7. the extended magneto hydrodynamics resistivity: etax
  8. 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
    flux2=np.zeros((n,nQ))
    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))
    st=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 [24]:
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]
            

Boundary conditions

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.

In [21]:
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))

Computation and outputs

Let's define some plotting parameters

In [22]:
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!

In [23]:
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")
DONE

© Pierre Gourdain, The University of Rochester, Charlie Seyler, Cornell University, 2020

In [ ]: