This is an automated email from the ASF dual-hosted git repository.
lahirujayathilake pushed a commit to branch agent-framewok-refactoring
in repository https://gitbox.apache.org/repos/asf/airavata.git
The following commit(s) were added to refs/heads/agent-framewok-refactoring by
this push:
new 93509da3a7 updated the notebook and docker files
93509da3a7 is described below
commit 93509da3a7c8b3ae4929f8c8e0dfceae9d8213e3
Author: lahiruj <[email protected]>
AuthorDate: Mon Mar 31 04:00:20 2025 -0400
updated the notebook and docker files
---
.../jupyterhub/user-container/MD/Dockerfile | 12 -
.../user-container/MD/build-container.sh | 4 +-
.../jupyterhub/user-container/MD/data/Allen/t.txt | 0
.../MD/data/Gkeyll/amitava_lectures_gpp2.ipynb | 1031 ++++++++++++++++++++
.../{plotE_z.ipynb => post_gkeyll_plotE_z.ipynb} | 38 +-
.../jupyterhub/user-container/MD/init.sh | 10 -
6 files changed, 1062 insertions(+), 33 deletions(-)
diff --git a/dev-tools/deployment/jupyterhub/user-container/MD/Dockerfile
b/dev-tools/deployment/jupyterhub/user-container/MD/Dockerfile
index 808cff06bd..a75a237739 100644
--- a/dev-tools/deployment/jupyterhub/user-container/MD/Dockerfile
+++ b/dev-tools/deployment/jupyterhub/user-container/MD/Dockerfile
@@ -2,22 +2,10 @@ FROM jupyter/base-notebook:latest
COPY data/ /tmp/default_data/
-COPY labconfig/jupyter_lab_config.py /jupyter_lab_config.py
-COPY labconfig/airavata_magics.py /airavata_magics.py
-COPY labconfig/__init__.py /__init__.py
-COPY labconfig/device_auth.py /device_auth.py
-COPY labconfig/bootstrap.sh /bootstrap.sh
-
COPY init.sh /usr/local/bin/init.sh
-# Create IPython startup directory
-RUN mkdir -p /home/jovyan/.ipython/profile_default/startup
-COPY labconfig/ipython_startup.py
/home/jovyan/.ipython/profile_default/startup/ipython_startup.py
-
USER root
RUN chmod +x /usr/local/bin/init.sh
-RUN chmod +x /bootstrap.sh
-RUN chown -R jovyan:users /home/jovyan/.ipython
USER $NB_USER
diff --git
a/dev-tools/deployment/jupyterhub/user-container/MD/build-container.sh
b/dev-tools/deployment/jupyterhub/user-container/MD/build-container.sh
index 4270dc6676..06e2e16d00 100755
--- a/dev-tools/deployment/jupyterhub/user-container/MD/build-container.sh
+++ b/dev-tools/deployment/jupyterhub/user-container/MD/build-container.sh
@@ -1,2 +1,2 @@
-docker build --platform linux/x86_64 -t cybershuttle/jupyterlab-base .
-docker push cybershuttle/jupyterlab-base
\ No newline at end of file
+docker build --platform linux/x86_64 -t cybershuttle/jupyterlab-base-sample .
+docker push cybershuttle/jupyterlab-base-sample
\ No newline at end of file
diff --git a/dev-tools/deployment/jupyterhub/user-container/MD/data/Allen/t.txt
b/dev-tools/deployment/jupyterhub/user-container/MD/data/Allen/t.txt
deleted file mode 100644
index e69de29bb2..0000000000
diff --git
a/dev-tools/deployment/jupyterhub/user-container/MD/data/Gkeyll/amitava_lectures_gpp2.ipynb
b/dev-tools/deployment/jupyterhub/user-container/MD/data/Gkeyll/amitava_lectures_gpp2.ipynb
new file mode 100644
index 0000000000..1517b2207d
--- /dev/null
+++
b/dev-tools/deployment/jupyterhub/user-container/MD/data/Gkeyll/amitava_lectures_gpp2.ipynb
@@ -0,0 +1,1031 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "67876c8f-bb64-4d65-a0e6-2eac8f5d81ee",
+ "metadata": {},
+ "source": [
+ "# Magnetic Reconnection Lectures Monday 03/31/25 and Wednesday
04/02/25\n",
+ "---\n",
+ "\n",
+ "**Goals:**\n",
+ "\n",
+ "\n",
+ "1. Reproduce the conservation arguments of the Sweet-Parker theory from
direct numerical simulations of nonlinear plasma fluid equations. \n",
+ "2. Examine what is meant by the concept of \"fast\" magnetic
reconnection in comparison to the original Sweet-Parker theory.\n",
+ "3. Determine what causes the reconnection to be \"fast\" by an Ohm's
Law analysis of our nonlinear simulations. \n",
+ "\n",
+ "---\n",
+ "\n",
+ "### Part 1: Some preliminaries\n",
+ "\n",
+ "To begin, let us review the equations we will numerically integrate to
model magnetic reconnection. We anticipate needing physics \"beyond MHD\" due
to the arguments laid out in last week's class on the Sweet-Parker theory:
while resistivity is sufficient non-ideal physics to break field lines, the
scaling of the inflow velocity, and thus the rate of magnetic reconnection, is
still too slow to explain observations of magnetic reconnection in laboratory,
space, and astrophysical plas [...]
+ "\n",
+ "We first consider the Vlasov-Maxwell system of equations\n",
+ "\\begin{align}\n",
+ " & \\frac{\\partial f_s}{\\partial t} + \\nabla \\cdot \\left
(\\mathbf{v} f_s \\right ) + \\nabla_{\\mathbf{v}} \\cdot \\left (
\\frac{q_s}{m_s} \\left [\\mathbf{E} + \\mathbf{v} \\times \\mathbf{B} \\right]
f_s \\right ) = C[f_s] \\\\\n",
+ " & \\epsilon_0 \\mu_0 \\frac{\\partial \\mathbf{E}}{\\partial t} -
\\nabla \\times \\mathbf{B} = \\mu_0 \\mathbf{J}, \\qquad \\mathbf{J} = \\sum_s
q_s \\int \\mathbf{v} f_s \\thinspace d\\mathbf{v} \\\\\n",
+ " & \\nabla \\cdot \\mathbf{E} = \\frac{\\rho_c}{\\epsilon_0}, \\qquad
\\rho_c = \\sum_s q_s \\int f_s \\thinspace d\\mathbf{v} \\\\\n",
+ " & \\frac{\\partial \\mathbf{B}}{\\partial t} + \\nabla \\times
\\mathbf{E} = 0, \\qquad \\nabla \\cdot \\mathbf{B} = 0 \n",
+ "\\end{align}\n",
+ "\n",
+ "To derive fluid equations from the Vlasov-Maxwell system, we take
(mass-weighted) velocity moments of the evolution equation for the distribution
function. The simplest sets of fluid equations utilize the zeroth, first, and
either scalar second or tensor second moment of the equation. The zeroth moment
gives the familiar conservation of mass equation: \n",
+ "\\begin{align}\n",
+ " \\frac{\\partial \\rho_s}{\\partial t} + \\nabla \\cdot \\left (
\\rho_s \\mathbf{u}_s \\right ) = 0 \n",
+ "\\end{align}\n",
+ "The first moment gives the familiar conservation of momentum equation \n",
+ "\\begin{align}\n",
+ " \\frac{\\partial \\rho_s \\mathbf{u}_s}{\\partial t} + \\nabla \\cdot
\\left ( \\rho_s \\mathbf{u}_s \\mathbf{u}_s + \\mathbf{P}_s \\right ) =
\\frac{q_s}{m_s} \\left ( \\rho_s \\mathbf{E} + \\rho_s \\mathbf{u}_s \\times
\\mathbf{B} \\right ) + \\mathbf{R}\n",
+ "\\end{align}\n",
+ "where $\\mathbf{P}_s$ is the pressure tensor of species $s$ and
$\\mathbf{R}$ is an as-yet-unspecified inter-species friction term. For at
least one system of equations we will solve today, we will take the
inter-species friction to be (for an electron-singly ionized ion plasma):\n",
+ "\\begin{align}\n",
+ " \\mathbf{R}_{ei} & = \\frac{\\rho_e}{\\tau_{ei}} \\left (\\mathbf{u}_e
- \\mathbf{u}_i \\right ), \\\\\n",
+ " \\tau_{ei} & = \\frac{6 \\sqrt{2} \\pi^{\\frac{3}{2}} \\epsilon_0^2 m_i
\\sqrt{m_e} T_e^{\\frac{3}{2}}}{\\ln \\Lambda e^4 \\rho_i}. \n",
+ "\\end{align}\n",
+ "If we assume $n_i \\approx n_e$, we can see how we can rearrange this
inter-species friction term into a more familiar form\n",
+ "\\begin{align}\n",
+ " \\mathbf{R}_{ei} & = \\frac{\\ln \\Lambda e^3 \\sqrt{m_e} n}{6
\\sqrt{2} \\pi^{\\frac{3}{2}} \\epsilon_0^2 T_e^{\\frac{3}{2}}} \\left
(en\\mathbf{u}_e - en\\mathbf{u}_i \\right ) \\\\\n",
+ " & = \\frac{\\ln \\Lambda}{6 \\sqrt{2} \\pi^{\\frac{3}{2}}}
\\frac{\\omega_{pe}}{e n \\lambda_{De}^3} \\left (en\\mathbf{u}_e -
en\\mathbf{u}_i \\right ) \\\\\n",
+ " & = -\\eta \\mathbf{J}\n",
+ "\\end{align}\n",
+ "**We will not assume $n_i \\approx n_e$ going forward; we have just
chosen to illustrate that this form of inter-species friction is a
generalization of the resistivity you have been working with up until this
point, where now the current density is given by the difference in flows
between the electrons and singly-charged ions**\n",
+ "\n",
+ "Now we consider the tensor second moment of the Vlasov equation,
switching to Einstein summation notation for clarity on the coupling to the
rank-3 heat-flux tensor and neglecting the contribution from the collision
operator: \n",
+ "\\begin{align}\n",
+ "\\frac{\\partial \\mathcal{P}_{ij_s}}{\\partial t} + \\frac{\\partial
\\mathcal{Q}_{ijk_s}}{\\partial x_k} = \\frac{q_s}{m_s} \\rho_s u_{s[i} E_{j]}
+ \\frac{q_s}{m_s} \\epsilon_{[ikl} \\mathcal{P}_{kj_s]} B_l\n",
+ "\\end{align}\n",
+ "where the square brackets around indices represent the minimal sum over
permutations of free indices needed to yield completely symmetric tensors. So
for example, $u_{s[i} E_{j]} = u_{i_s} E_j + u_{j_s} E_i$. This tensor second
moment $\\mathcal{P}_{ij_s} = P_{ij_s} + \\rho_s u_{i_s} u_{j_s} =
\\mathbf{P}_s + \\rho_s \\mathbf{u}_s \\mathbf{u}_s$ is often referred to as
the stress tensor and includes both the pressure tensor and Reynolds stress.
The rank-3 heat-flux tensor can be wri [...]
+ "\\begin{align}\n",
+ " \\mathcal{Q}_{ijk_s} = Q_{ijk_s} + u_{s[i} \\mathcal{P}_{jk_s]} - 2
\\rho_s u_{i_s} u_{j_s} u_{k_s}\n",
+ "\\end{align}\n",
+ "where\n",
+ "\\begin{align}\n",
+ " Q_{ijk_s} = \\int (v_i - u_{i_s}) (v_j - u_{j_s}) (v_k - u_{k_s}) f_s
\\thinspace d\\mathbf{v}\n",
+ "\\end{align}\n",
+ "is the plasma rest-frame heat flux. We now see what is commonly referred
to as the \"closure\" problem as we have 10 equations, 1 from the mass density
evolution, 3 from the momentum density evolution, and 6 from the stress tensor
evolution, but 20 unknowns thanks to the 10 unknown components for the plasma
rest-frame heat flux. To close the system, we must express the plasma
rest-frame heat flux in terms of lower moments. A straightforward equation we
will use expresses \n",
+ "\\begin{align}\n",
+ " \\frac{\\partial Q_{ijk_s}}{\\partial x_k} & = v_{th_s} k_0 \\left
(P_{ij_s} - p_s \\right ) \\\\ \n",
+ " p_s & = \\frac{1}{3} P_{ii}, \\qquad v_{th_s} =
\\sqrt{\\frac{p_s}{\\rho_s}}.\n",
+ "\\end{align}\n",
+ "In other words, the plasma rest-frame heat flux works to relax the
pressure tensor back to isotropy and a scalar pressure at a specified length
scale $k_0$. \n",
+ "\n",
+ "Based on the form of the tensor second moment equation, or directly from
the Vlasov equation **neglecting collisions**, what is the scalar second moment
equation for the evolution of the total particle energy?\n",
+ "1. $\\frac{\\partial \\mathcal{E}_s}{\\partial t} + \\nabla \\cdot
\\left ( \\left [\\mathcal{E}_s + p_s\\right ] \\mathbf{u}_s + \\mathbf{q}_s
\\right ) = \\frac{q_s}{m_s} \\rho_s \\mathbf{u}_s \\cdot \\mathbf{E} $\n",
+ "2. $\\frac{\\partial \\mathcal{E}_s}{\\partial t} + \\nabla \\cdot
\\left ( \\left [\\mathcal{E}_s + p_s\\right ] \\mathbf{u}_s + \\mathbf{q}_s
\\right ) = \\frac{q_s}{m_s} \\left (\\rho_s \\mathbf{u}_s \\cdot \\mathbf{E} +
\\rho_s |\\mathbf{u}_s \\times \\mathbf{B}|^2 \\right )$\n",
+ "3. $\\frac{\\partial \\mathcal{E}_s}{\\partial t} + \\nabla \\cdot
\\left ( \\left [\\mathcal{E}_s + p_s\\right ] \\mathbf{u}_s + \\mathbf{q}_s
\\right ) = 0$\n",
+ "4. $\\frac{\\partial \\mathcal{E}_s}{\\partial t} + \\nabla \\cdot
\\left ( \\left [\\mathcal{E}_s + p_s\\right ] \\mathbf{u}_s + \\mathbf{q}_s
\\right ) = \\mathbf{J} \\cdot \\left(\\mathbf{u}_s \\times \\mathbf{B}\\right
) $\n",
+ "\n",
+ "[Answer Form
1.](https://docs.google.com/forms/d/e/1FAIpQLSetFle7blYJSFAuEqt2WItk0KZ_q5v9q1-OPfkOrbYYLTDKCA/viewform?usp=dialog)\n",
+ "\n",
+ "To close the fluid equations at this conservation of energy equation,
assuming an ideal fluid, we can set $\\mathbf{q}_s = 0$ and obtain the scalar
pressure via the equation of state, $\\mathcal{E}_s = \\frac{p_s}{\\Gamma - 1}
+ \\frac{1}{2} \\rho_s |\\mathbf{u}_s|^2$, where $\\Gamma$ is the adiabatic
index of the fluid. "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "f232b540-9a34-4b3e-9f55-06d860e87d89",
+ "metadata": {},
+ "source": [
+ "### Some terminology and helper functions\n",
+ "\n",
+ "The coupled set of fluid equations: either the 5 equations from our
velocity moments through a scalar second moment, or the 10 equations through
the tensor second velocity moment, **we refer to as the 5m and 10m, or 5 moment
and 10 moment, multi-fluid equations.** We can solve these equations for any
number of species (even mixing some fluids being 5 moment while other fluids
are 10 moment) and couple to Maxwell's equations via the current density, which
is simply a sum over species [...]
+ "\n",
+ "We can now proceed with the necessary preliminaries, including importing
our needed Python libraries such as matplotlib and numpy, and defining some
helper functions which will be useful for visualizing and analyzing the data. "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "9fcdfc96-f11e-4820-b385-a888dfe1c42e",
+ "metadata": {},
+ "source": [
+ "#### Initialize, Authenticate, and Configure Remote Execution
Environment\n",
+ "\n",
+ "##### Start a Cybershuttle Remote Executor: This will submit a remote job
to Anvil\n",
+ "- **Cluster** - Anvil & Jetstream\n",
+ "- **Community Allocation** - Gkeyll\n",
+ "- **Requested Resources** \n",
+ " - Anvil:4 CPUs, 4GB Memory, 60 Minutes, Shared Queue\n",
+ " - Jetstream:1 CPUs, 60 Minutes, Cloud Queue\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "515df1fa-13dd-43f7-83c8-1fd118f4865b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "%pip install -q --force-reinstall airavata-python-sdk[notebook]\n",
+ "import airavata_jupyter_magic\n",
+ "\n",
+ "%authenticate\n",
+ "\n",
+ "# PLEASE NOTE: At a given time, ONLY run a single job in one of the
available clusters: Anvil OR Jetstream. \n",
+ "#Anvil is the default cluster and if you need to run on Jetstream;
comment the Anvil and uncomment the Jetstream\n",
+ "\n",
+ "#Running on Anvil remote cluster\n",
+ "%request_runtime test_cpu --cluster=Anvil --cpus=4 --memory=4096
--walltime=60 --queue=shared --group=Gkeyll\n",
+ "%switch_runtime test_cpu\n",
+ "\n",
+ "#Running on Jetstream cluster\n",
+ "#%request_runtime test_cpu --cluster=JS --cpus=1 --walltime=60
--queue=cloud --group=Gkeyll\n",
+ "#%switch_runtime test_cpu"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "f189ff7b-db4f-47fe-affa-66f1182f25c1",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "#Defininf the data paths in remote clusters\n",
+ "%ls -l /data\n",
+ "\n",
+ "#Anvil and Jetstream data path\n",
+ "data_path = \"/data/amitava-class\"\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "f34973fc-90bd-45d5-9ce6-6ff24af0a8b8",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import matplotlib\n",
+ "import matplotlib.pyplot as plt\n",
+ "import postgkyl as pg\n",
+ "import numpy as np\n",
+ "\n",
+ "#Note: This cell may execute longer if the job is not yet active in the
remote cluster. Please pause till the status =READY message appears."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "5e4c1ef9-7c46-4449-a977-a3029e18d3d6",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Helper function for computing vector potential\n",
+ "def calc_psi2d(fx, fy, dx=1, dy=1):#solenoidal flows\n",
+ "\t'''\n",
+ "\tCalcualte psi by integrating dpsi = -fy*dx + fx*dy, psi[0,0]=0.\n",
+ "\tNotes: \n",
+ "\t\t1. (fx=dpsi/dy,fy=-dpsi/dx) is called Hamiltonian gradient of psi,
and\tcontours of psi give vector field (fx, fy);\n",
+ "\t\t2. div(f)=0\n",
+ "\t'''\n",
+ "\tny,nx=fx.shape\n",
+ "\tpsi=np.zeros((ny,nx))\n",
+ "\tfor jx in range(1,nx):\n",
+ "\t\tpsi[0,jx]=psi[0,jx-1]-fy[0,jx]*dx\n",
+ "\tfor jy in range(1,ny):\n",
+ "\t\tpsi[jy,:]=psi[jy-1,:]+fx[jy,:]*dy\n",
+ "\t# since f = rot(A) gives extra restraints on f (e.g., div(f)=0)\n",
+ "\t# it makes sense that information provided by fy[1:nx,:] is useless
here\n",
+ "\treturn psi\n",
+ "print(f\"{data_path}\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "873222c7-359a-4184-9862-c4e448fa87f6",
+ "metadata": {},
+ "source": [
+ "### Defining our units\n",
+ "\n",
+ "When running simulations, it is often convenient to define a set of
normalized units in which to perform the simulation. Many physical constants,
such as the permittivity of free space, $\\epsilon_0 = 8.85 \\times 10^{-12}
F/m$, or the actual electron and proton mass, $m_e = 9.11\\times 10^{-31} kg,
m_p = 1.67\\times 10^{-27}kg$ respectively, could lead to enhanced truncation
error due to finite precision arithmetic if we carried around the physical
values of these constants in our [...]
+ "\n",
+ "Below I define the units of the electron-ion multi-fluid magnetic
reconnection simulations we have performed. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4d5bbc3a-bb3d-4743-bd2d-60d9a7f3ab7d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Physical constants and derived parameters\n",
+ "gas_gamma = 5.0 / 3.0 # Adiabatic index.\n",
+ "epsilon0 = 1.0 # Permittivity of free space.\n",
+ "mu0 = 1.0 # Permeability of free space.\n",
+ "light_speed = 1.0/np.sqrt(epsilon0*mu0) # Speed of light. \n",
+ "mass_ion = 100.0 # Ion mass.\n",
+ "charge_ion = 1.0 # Ion charge.\n",
+ "mass_elc = 1.0 # Electron mass.\n",
+ "charge_elc = -1.0 # Electron charge.\n",
+ "Ti_over_Te = 5.0 # Ion temperature / electron temperature.\n",
+ "n0 = 1.0 # Reference number density.\n",
+ "nb_over_n0 = 0.2 # Background number density / reference number
density.\n",
+ "wpi = np.sqrt(charge_ion**2 * n0 / (epsilon0 * mass_ion)) # Ion plasma
frequency. \n",
+ "wpe = np.sqrt(charge_ion**2 * n0 / (epsilon0 * mass_elc)) # Electron
plasma frequency. \n",
+ "di = light_speed/wpi # Ion inertial length. \n",
+ "de = light_speed/wpe # Electron inertial length. \n",
+ "B0 = 0.5 # Reference magnetic field strength. Sets ratio of vAe/c. \n",
+ "beta = 1.0 # Plasma beta.\n",
+ "omega_ci = charge_ion * B0 / mass_ion # Ion cyclotron frequency. \n",
+ "vA0 = B0/np.sqrt(mu0*n0*mass_ion) # Reference Alfven speed. \n",
+ "vA_up = B0/np.sqrt(mu0*nb_over_n0*mass_ion) # Upstream reference Alfven
speed from background density. \n",
+ "Lambda_ee = 2.0e5 # Plasma parameter for computing collision time. \n",
+ "Ti_frac = Ti_over_Te / (1.0 + Ti_over_Te) # Fraction of total temperature
from ions.\n",
+ "Te_frac = 1.0 / (1.0 + Ti_over_Te) # Fraction of total temperature from
electrons.\n",
+ "T_tot = beta * (B0 * B0) / 2.0 / n0 # Total temperature.\n",
+ "T_elc_ref = T_tot*Te_frac # Reference electron temperature. "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "dfb5a86e-38ae-430f-a176-9630d8a03ebc",
+ "metadata": {},
+ "source": [
+ "After examination of these units, my normalized time, length, and
velocity scales are?\n",
+ "1. $t = \\omega_{pe}^{-1} \\tilde{t} $, $\\mathbf{x} = d_e
\\tilde{\\mathbf{x}}$, $\\mathbf{v} = c \\tilde{\\mathbf{v}}$\n",
+ "2. $t = \\omega_{pi}^{-1} \\tilde{t} $, $\\mathbf{x} = d_i
\\tilde{\\mathbf{x}}$, $\\mathbf{v} = c \\tilde{\\mathbf{v}}$\n",
+ "3. $t = \\omega_{pe}^{-1} \\tilde{t} $, $\\mathbf{x} = \\lambda_{De}
\\tilde{\\mathbf{x}}$, $\\mathbf{v} = v_{th_e} \\tilde{\\mathbf{v}}$\n",
+ "4. $t = \\Omega_{ci}^{-1} \\tilde{t} $, $\\mathbf{x} = \\rho_i
\\tilde{\\mathbf{x}}$, $\\mathbf{v} = v_{th_i} \\tilde{\\mathbf{v}}$\n",
+ "\n",
+ "where $\\omega_{ps} = \\sqrt{e^2 n_s/\\epsilon_0 m_s}$ is the plasma
frequency of species $s$, $d_s = c/\\omega_{ps}$ is the inertial length of
species $s$, $\\lambda_{De} = v_{th_e}/\\omega_{pe}$ is the Debye length,
$\\Omega_{ci} = e B_0/m_i$ is the ion cyclotron frequency, and $\\rho_i =
\\sqrt{2} v_{th_i}/\\Omega_{ci}$ is the ion gyroradius. \n",
+ "\n",
+ "If we include resistivity in our 5 moment calculations, we need to
specify the plasma parameter, $\\Lambda \\sim n \\lambda_{De}^3$, to compute
the collision time $\\tau_{ei}$. In these normalized units, the collision time
scales like\n",
+ "1. $\\tau_{ei} \\sim \\Lambda \\tilde{t}$\n",
+ "2. $\\tau_{ei} \\sim \\ln\\Lambda \\tilde{t}$\n",
+ "3. $\\tau_{ei} \\sim \\frac{\\Lambda}{\\ln\\Lambda} \\tilde{t}$\n",
+ "4. $\\tau_{ei} \\sim \\frac{\\ln\\Lambda}{\\Lambda} \\tilde{t}$\n",
+ "\n",
+ "[Answer Form
2](https://docs.google.com/forms/d/e/1FAIpQLSflCoqOnP96nbp_jnW4KR6HOe6f4atx7sFtP2FUMG782WxErQ/viewform?usp=sharing)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "c9e43d49-dacf-4f3c-acef-1888c628643d",
+ "metadata": {},
+ "source": [
+ "### $\\texttt{Gkeyll}$ and $\\texttt{postgkyl}$\n",
+ "\n",
+ "The code we will be using to model magnetic reconnection is
[$\\texttt{Gkeyll}$](https://gkeyll.readthedocs.io/), a general purpose
simulation framework for a variety of fluid and plasma systems. You can
download and install $\\texttt{Gkeyll}$ yourself by following the installation
instructions on our [Github repo](https://github.com/ammarhakim/gkylzero). \n",
+ "\n",
+ "To read the data, we will utilize the post-processing suite we have
developed alongside $\\texttt{Gkeyll}$,
[$\\texttt{postgkyl}$](https://github.com/ammarhakim/postgkyl), which you can
also download and install via the instructions on Github. The cluster we will
be utilizing for analyzing the results of our simulations already has
installations of $\\texttt{Gkeyll}$ and $\\texttt{postgkyl}$; we already
imported postgkyl in this Jupyter Notebook, so if we did not have
$\\texttt{post [...]
+ "\n",
+ "The output of $\\texttt{Gkeyll}$ simulations can be manipulated in one of
two ways: through the GData class, which retains useful metadata from the
simulation to subsequent operations, or by directly fetching the raw values and
grid and storing them in Numpy arrays for our subsequent manipulations. We will
utilize both means of reading the data as we analyze our simulation results. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "87f022c1-770a-416b-be68-76a49f3ee084",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def get_Gdata(pre, frame):\n",
+ " data_elc = pg.data.GData(\"%s-elc_%d.gkyl\" % (pre, frame))\n",
+ " data_ion = pg.data.GData(\"%s-ion_%d.gkyl\" % (pre, frame))\n",
+ " data_field = pg.data.GData(\"%s-field_%d.gkyl\" % (pre, frame))\n",
+ "\n",
+ " return data_elc, data_ion, data_field"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "6ea76d1f-d510-4b18-a166-554acc48f6d0",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def read_data(data_elc, data_ion, data_field):\n",
+ " fluid_elc = data_elc.get_values()\n",
+ " fluid_ion = data_ion.get_values() \n",
+ " field = data_field.get_values() \n",
+ "\n",
+ " # Same grid for all electrons, ions, and EM fields\n",
+ " coords = data_elc.get_grid()\n",
+ " # Center the grid values\n",
+ " for d in range(len(coords)):\n",
+ " coords[d] = 0.5*(coords[d][:-1] + coords[d][1:])\n",
+ "\n",
+ " return coords, fluid_elc, fluid_ion, field"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "8eb6eedf-9e1e-4b3a-a792-d9b44bcb6714",
+ "metadata": {},
+ "source": [
+ "### A look at the initial conditions\n",
+ "\n",
+ "In the beginning of the tearing mode calculation we started last week, we
considered an equilibrium magnetic field profile of the form\n",
+ "\\begin{align}\n",
+ "\\mathbf{B}(x,y) = B_0 \\tanh \\left ( \\frac{y}{\\lambda} \\right )
\\hat{\\mathbf{x}} \n",
+ "\\end{align}\n",
+ "Such a magnetic field necessarily has a corresponding current density
profile of the form\n",
+ "\\begin{align}\n",
+ " \\mathbf{J}(x,y) = -\\frac{B_0}{\\lambda} \\textrm{sech}^2 \\left (
\\frac{y}{\\lambda} \\right ) \\hat{\\mathbf{z}}\n",
+ "\\end{align}\n",
+ "Is this a force-free equilibrium? Can we re-write this equilibrium as
$\\nabla \\times \\mathbf{B} = \\alpha \\mathbf{B}$ for some $\\alpha$?\n",
+ "1. Yes\n",
+ "2. No\n",
+ "\n",
+ "If this equilibrium is not force-free, what pressure profile do we need
to satisfy force-balance? \n",
+ "1. No pressure is needed; this equilibrium is force-free\n",
+ "2. $p \\propto \\textrm{sech}^2 \\left ( \\frac{y}{\\lambda} \\right
)$\n",
+ "3. $p \\propto \\tanh \\left ( \\frac{y}{\\lambda} \\right )
\\textrm{sech}^2 \\left ( \\frac{y}{\\lambda} \\right )$\n",
+ "4. $p \\propto \\tanh \\left ( \\frac{y}{\\lambda} \\right )$\n",
+ "\n",
+ "[Answer Form
3](https://docs.google.com/forms/d/e/1FAIpQLScw_ZUM63GHiY4jdReTyHAmrbvGpClxp-R5vn1knsUDhlEZ1w/viewform?usp=dialog)\n",
+ "\n",
+ "Let's visualize our equilibrium from our simulations. Since the
equilibrium only requires the lowest moments, the equilibrium is identical in
both 5-moment and 10-moment. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "aa929847-7e22-468a-97d4-1140625650b5",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "data_elc, data_ion, data_field =
get_Gdata(f\"{data_path}/resistive_gem_mass100/rt_5m_resistive_gem\", 0)\n",
+ "r_coords, r_fluid_elc, r_fluid_ion, r_field = read_data(data_elc,
data_ion, data_field)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "06bb4782-e485-44a4-9727-1fc001d6f085",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Compute out-of-plane current density\n",
+ "Nx = r_fluid_ion.shape[0]\n",
+ "Ny = r_fluid_ion.shape[1]\n",
+ "Jz = np.zeros((Nx, Ny))\n",
+ "Jz = charge_ion/mass_ion*r_fluid_ion[...,3] +
charge_elc/mass_elc*r_fluid_elc[...,3]\n",
+ "# Compute contours of the magnetic field.\n",
+ "Bx = r_field[...,3]\n",
+ "By = r_field[...,4]\n",
+ "dx = r_coords[0][1] - r_coords[0][0]\n",
+ "dy = r_coords[1][1] - r_coords[1][0]\n",
+ "psi = np.zeros((Ny, Nx))\n",
+ "psi = calc_psi2d(Bx.transpose(),By.transpose(), dx, dy)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d5a7ab12-9b7c-4cf9-b078-6df69c6d665b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.figure()\n",
+ "# Add shading='gouraud' for further interpolation and smoother plots\n",
+ "plt.pcolormesh(r_coords[0]/di, r_coords[1]/di, Bx.transpose(), vmax=-0.5,
vmin=0.5, cmap='inferno')#, shading='gouraud')\n",
+ "plt.colorbar(format='%.2f', ticks=np.linspace(-0.5, 0.5, 5),
fraction=0.046, pad=0.04)\n",
+ "plt.xlabel(r'$x/d_i$')\n",
+ "plt.ylabel(r'$y/d_i$')\n",
+ "plt.title(r'$B_x$')\n",
+ "plt.setp(plt.gca(), aspect=1.0)\n",
+ "\n",
+ "vmax = np.max(np.abs(Jz))\n",
+ "vmin = -vmax\n",
+ "plt.figure()\n",
+ "# Add shading='gouraud' for further interpolation and smoother plots\n",
+ "plt.pcolormesh(r_coords[0]/di, r_coords[1]/di, Jz.transpose(), vmax=vmax,
vmin=vmin, cmap='seismic')#, shading='gouraud')\n",
+ "plt.colorbar(format='%.2f', ticks=np.linspace(vmin, vmax, 5),
fraction=0.046, pad=0.04)\n",
+ "plt.contour(r_coords[0]/di, r_coords[1]/di, psi, 7, colors=\"k\",
linestyles=\"solid\")\n",
+ "plt.xlabel(r'$x/d_i$')\n",
+ "plt.ylabel(r'$y/d_i$')\n",
+ "plt.title(r'$J_z$')\n",
+ "plt.setp(plt.gca(), aspect=1.0)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "8a87217d-ae4a-4fd3-9824-eea12ea8c2e7",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Plot Bx and Jz at x=0\n",
+ "x0 = int(Nx/2)\n",
+ "plt.figure()\n",
+ "plt.plot(r_coords[1]/di, Bx[x0,...], 'k', label=r'$B_x$')\n",
+ "plt.plot(r_coords[1]/di, Jz[x0,...], 'b', label=r'$J_z$')\n",
+ "plt.xlabel(r'$y/d_i$')\n",
+ "plt.xlim(-6.0, 6.0)\n",
+ "plt.ylim(-0.5, 0.5)\n",
+ "plt.legend(loc='best')\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "1fbdc4b1-98d1-4cc0-8108-9f091dd33e89",
+ "metadata": {},
+ "source": [
+ "### A look later in time\n",
+ "\n",
+ "Let's look at the simulation later in time to see the fully developed
reconnection. You can feel free to switch between visualizing the 10 moment and
5 moment data, or even adding your own calls to the read and plot routines to
visualize the development in time. \n",
+ "\n",
+ "**We have written out the data every $\\Delta t_{IO} = 0.5
\\Omega_{ci}^{-1}$ to time $t_{end} = 30 \\Omega_{ci}^{-1}$ for a total of 60
frames**."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "29f0c229-8892-4938-a4d3-52bc4cdcbed4",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "#data_elc, data_ion, data_field =
get_Gdata(f\"{data_path}/10m_gem_mass100/rt_10m_gem\", 60)\n",
+ "data_elc, data_ion, data_field =
get_Gdata(f\"{data_path}/resistive_gem_mass100/rt_5m_resistive_gem\", 60)\n",
+ "r_coords, r_fluid_elc, r_fluid_ion, r_field = read_data(data_elc,
data_ion, data_field)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "322cae84-ed8c-45e7-aca2-359e4f13e1a3",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Compute out-of-plane current density\n",
+ "Nx = r_fluid_ion.shape[0]\n",
+ "Ny = r_fluid_ion.shape[1]\n",
+ "Jz = np.zeros((Nx, Ny))\n",
+ "Jz = charge_ion/mass_ion*r_fluid_ion[...,3] +
charge_elc/mass_elc*r_fluid_elc[...,3]\n",
+ "# Compute contours of the magnetic field.\n",
+ "Bx = r_field[...,3]\n",
+ "By = r_field[...,4]\n",
+ "dx = r_coords[0][1] - r_coords[0][0]\n",
+ "dy = r_coords[1][1] - r_coords[1][0]\n",
+ "psi = np.zeros((Ny, Nx))\n",
+ "psi = calc_psi2d(Bx.transpose(),By.transpose(), dx, dy)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "26853cb3-bda2-45b4-aeba-be82d989d68b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "vmax = np.max(np.abs(Jz))\n",
+ "vmin = -vmax\n",
+ "plt.figure()\n",
+ "# Add shading='gouraud' for further interpolation and smoother plots\n",
+ "plt.pcolormesh(r_coords[0]/di, r_coords[1]/di, Jz.transpose(), vmax=vmax,
vmin=vmin, cmap='seismic')#, shading='gouraud')\n",
+ "plt.colorbar(format='%.1f', ticks=np.linspace(vmin, vmax, 7),
fraction=0.046, pad=0.04)\n",
+ "plt.contour(r_coords[0]/di, r_coords[1]/di, psi, 7, colors=\"k\",
linestyles=\"solid\")\n",
+ "plt.xlabel(r'$x/d_i$')\n",
+ "plt.ylabel(r'$y/d_i$')\n",
+ "plt.title(r'$J_z$')\n",
+ "plt.setp(plt.gca(), aspect=1.0)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "4be5854f-2219-448e-9abc-91332ab066ea",
+ "metadata": {},
+ "source": [
+ "### Determining the scaling of the in-flows and out-flows: how good is
the Sweet-Parker argument?\n",
+ "\n",
+ "We now will address the first of our goals today: reconstructing the
conservation arguments which Sweet and Parker utilized originally to argue for
magnetic reconnection as a possible energy conversion process in observed
phenomena such as solar flares. First let us confirm the simulation is behaving
as we expect; is the magnetic field vanishing at the X-point? Can we
intuitively see that some energy conversion process must be occurring?"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "ae8bcb38-bb78-4d3b-bd41-35fe0e0afb67",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Compute |B|^2 and |B|\n",
+ "_, magB_sq = pg.tools.mag_sq(data_field, '3:6')\n",
+ "magB = np.sqrt(magB_sq[...,0])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "c4beba44-1ad9-461a-ab29-b7fd4384f8f6",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "vmax = np.max(magB)\n",
+ "vmin = 0.0\n",
+ "plt.figure()\n",
+ "# Add shading='gouraud' for further interpolation and smoother plots\n",
+ "plt.pcolormesh(r_coords[0]/di, r_coords[1]/di, magB.transpose(),
vmax=vmax, vmin=vmin, cmap='inferno')#, shading='gouraud')\n",
+ "plt.colorbar(format='%.1f', ticks=np.linspace(vmin, vmax, 7),
fraction=0.046, pad=0.04)\n",
+ "plt.contour(r_coords[0]/di, r_coords[1]/di, psi, 7, colors=\"k\",
linestyles=\"solid\")\n",
+ "plt.xlabel(r'$x/d_i$')\n",
+ "plt.ylabel(r'$y/d_i$')\n",
+ "plt.title(r'$|\\mathbf{B}|$')\n",
+ "plt.setp(plt.gca(), aspect=1.0)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "8df6bdbc-4819-49dd-b2ff-7e8564c8f237",
+ "metadata": {},
+ "source": [
+ "Since the Alfvén speed is a useful normalization when examining the flows
produced by the reconnection, let's visualize the local Alfvén speed now from
this local magnetic field magnitude. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "cc7fc14b-3ad3-405e-bd35-203518f3fb5b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Compute vA = |B|/sqrt(mu_0 rho_i)\n",
+ "_, rho_ion = pg.tools.get_density(data_ion)\n",
+ "vA = magB/np.sqrt(mu0*rho_ion[...,0])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "325e2768-6d1d-422c-abe6-7ba8d53e7984",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "vmax = np.max(vA)\n",
+ "vmin = 0.0\n",
+ "plt.figure()\n",
+ "# Add shading='gouraud' for further interpolation and smoother plots\n",
+ "plt.pcolormesh(r_coords[0]/di, r_coords[1]/di, vA.transpose(), vmax=vmax,
vmin=vmin, cmap='inferno')#, shading='gouraud')\n",
+ "plt.colorbar(format='%.1f', ticks=np.linspace(vmin, vmax, 4),
fraction=0.046, pad=0.04)\n",
+ "plt.contour(r_coords[0]/di, r_coords[1]/di, psi, 7, colors=\"k\",
linestyles=\"solid\")\n",
+ "plt.xlabel(r'$x/d_i$')\n",
+ "plt.ylabel(r'$y/d_i$')\n",
+ "plt.title(r'$v_A$')\n",
+ "plt.setp(plt.gca(), aspect=1.0)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "b64f5607-0c7d-4c7b-97ac-48183f1ae832",
+ "metadata": {},
+ "source": [
+ "We have a couple different options for normalizing our in-flow and
out-flow velocities. Let's visualize both! "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "e9d4b23b-6f93-4740-bdac-ef465a5b4da4",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "coords, ux = pg.tools.get_vx(data_ion)\n",
+ "ux_norm0 = ux[...,0]/vA_up # Normalized x flow to upstream Alfven
speed\n",
+ "ux_norm = ux[...,0]/vA # Normalized x flow to local Alfven speed"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "150abc63-fa9a-49ad-a036-7f49077ba09b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "vmax = np.max(np.abs(ux_norm0))\n",
+ "vmin = -vmax\n",
+ "plt.figure()\n",
+ "# Add shading='gouraud' for further interpolation and smoother plots\n",
+ "plt.pcolormesh(r_coords[0]/di, r_coords[1]/di, ux_norm.transpose(),
vmax=vmax, vmin=vmin, cmap='inferno')#, shading='gouraud')\n",
+ "plt.colorbar(format='%.1f', ticks=np.linspace(vmin, vmax, 5),
fraction=0.046, pad=0.04)\n",
+ "plt.contour(r_coords[0]/di, r_coords[1]/di, psi, 7, colors=\"k\",
linestyles=\"solid\")\n",
+ "plt.xlabel(r'$x/d_i$')\n",
+ "plt.ylabel(r'$y/d_i$')\n",
+ "plt.title(r'$u_x/v_A$')\n",
+ "plt.setp(plt.gca(), aspect=1.0)\n",
+ "\n",
+ "plt.figure()\n",
+ "# Add shading='gouraud' for further interpolation and smoother plots\n",
+ "plt.pcolormesh(r_coords[0]/di, r_coords[1]/di, ux_norm0.transpose(),
vmax=vmax, vmin=vmin, cmap='inferno')#, shading='gouraud')\n",
+ "plt.colorbar(format='%.1f', ticks=np.linspace(vmin, vmax, 5),
fraction=0.046, pad=0.04)\n",
+ "plt.contour(r_coords[0]/di, r_coords[1]/di, psi, 7, colors=\"k\",
linestyles=\"solid\")\n",
+ "plt.xlabel(r'$x/d_i$')\n",
+ "plt.ylabel(r'$y/d_i$')\n",
+ "plt.title(r'$u_x/v_{A_{up}}$')\n",
+ "plt.setp(plt.gca(), aspect=1.0)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "6d736802-b3be-423c-90b6-72353dcef081",
+ "metadata": {},
+ "source": [
+ "And let's also take cuts of these quantities along $y=0$ to see the ion
out-flow, $u_x$, along the X-line. We can also plot the out-of-plane current
density, $J_z$, along $y=0$ to estimate the length of the reconnecting current
sheet, $L$. Is this what we expect? What should the outflow be based on Sweet
and Parker's conservation of energy argument?\n",
+ "1. $u_x \\sim V_{A_{x=y=0}} =
\\frac{|\\mathbf{B}|_{x=y=0}}{\\sqrt{\\mu_0 \\rho_{i_{x=y=0}}}}$\n",
+ "2. $u_x \\sim \\frac{|\\mathbf{B}|_{x=y=0}}{\\sqrt{\\mu_0
\\rho_{i_{up}}}}$\n",
+ "3. $u_x \\sim V_{A_{up}} = \\frac{B_0}{\\sqrt{\\mu_0
\\rho_{i_{up}}}}$\n",
+ "4. $u_x \\sim \\frac{B_0}{\\sqrt{\\mu_0 \\rho_{i_{x=y=0}}}}$\n",
+ "\n",
+ "Similarly, let's examine the in-flow velocity, $u_y$, and out-of-plane
current density, $J_z$, at $x=0$ to determine the width, $\\delta$, of the
current sheet. Is this also what we expect? What should the inflow be base on
Sweet and Parker's conservation of mass argument? \n",
+ "1. $u_y \\sim \\frac{L}{\\delta} u_x$\n",
+ "2. $u_y \\sim \\frac{\\delta}{L} u_x$\n",
+ "3. $u_y \\sim \\delta u_x $\n",
+ "4. $u_y \\sim L u_x$\n",
+ "\n",
+ "[Answer Form
4](https://docs.google.com/forms/d/e/1FAIpQLScSG_2YGwhf5hvFoluxxzPsxN13FXdHqlw0n5axdwx-xg9GQg/viewform?usp=dialog)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "99ebcfe8-2a30-414c-81d3-35e7dc46cab0",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Take a cut of ux/vA_up at y=0\n",
+ "ux_norm0_y0 = ux_norm0[..., int(Ny/2)]\n",
+ "Jz_y0 = Jz[..., int(Ny/2)]\n",
+ "plt.figure()\n",
+ "plt.plot(r_coords[0]/di, ux_norm0_y0, 'k')\n",
+ "plt.xlabel(r'$x/d_i$')\n",
+ "plt.ylabel(r'$u_x/v_{A_{up}}$')\n",
+ "plt.xlim(-12,12)\n",
+ "\n",
+ "# Plot Jz extent at y=0 to see length of current sheet\n",
+ "plt.figure()\n",
+ "plt.plot(r_coords[0]/di, Jz_y0, 'b')\n",
+ "plt.xlabel(r'$x/d_i$')\n",
+ "plt.ylabel(r'$J_z$')\n",
+ "plt.xlim(-5,5)\n",
+ "\n",
+ "# Fetch uy and Jz and take a cut through x = 0 to see width of current
sheet\n",
+ "coords, uy = pg.tools.get_vy(data_ion)\n",
+ "uy_x0 = uy[int(Nx/2), ..., 0]/vA_up\n",
+ "Jz_x0 = Jz[int(Nx/2), ...]\n",
+ "plt.figure()\n",
+ "plt.plot(r_coords[1]/di, uy_x0, 'k', label=r'$u_y/v_{A_{up}}$')\n",
+ "plt.plot(r_coords[1]/di, Jz_x0, 'b', label=r'$J_z$')\n",
+ "plt.xlabel(r'$y/d_i$')\n",
+ "plt.xlim(-6,6)\n",
+ "plt.legend(loc='best')\n",
+ "\n",
+ "plt.figure()\n",
+ "plt.plot(r_coords[1]/di, uy_x0, 'k', label=r'$u_y/v_{A_{up}}$')\n",
+ "plt.plot(r_coords[1]/di, Jz_x0, 'b', label=r'$J_z$')\n",
+ "plt.xlabel(r'$y/d_i$')\n",
+ "plt.xlim(-1,1)\n",
+ "plt.legend(loc='best')\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "7ee3a642-f35c-450f-9e35-221a63518b4e",
+ "metadata": {},
+ "source": [
+ "### So what does it even mean for the reconnection to be fast?\n",
+ "\n",
+ "We have now shown that our simulations obey the Sweet-Parker scalings.
Good! \\texttt{Gkeyll} obeys conservation of mass and energy! But wait a
minute, when we derived the Sweet-Parker scaling before we found the
reconnection rate before, it scaled with the resistivity, meaning the in-flow
velocity should be related to the Alfvén velocity by the resistivity. We know
that, \n",
+ "\\begin{align}\n",
+ "u_y \\sim u_x S^{-\\frac{1}{2}}, \\quad S = \\frac{L V_A}{\\eta}\n",
+ "\\end{align}\n",
+ "Let's estimate this from our resistive 5 moment calculation. Given the
choice of plasma parameter, $\\Lambda$, and $L \\sim 10 d_i$, we should get a
Lundquist number of $\\sim 10^5$ (try plugging in all the values for the code
units; since the Lundquist number is a normalized quantity you can just use the
code unit values for everything, $L, V_A, \\eta$, and safely estimate the
Lundquist number). \n",
+ "\n",
+ "Well.. that is much slower than we are finding right? So what is the
actual reconnection rate? We can compute this by either evaluating the
reconnecting electric field $E_z$ at the X-point, or computing the rate of
change of the vector potential. Either way, we see that in fact the rate is
actually $\\sim 0.1 V_A$, **much faster** than what Sweet and Parker estimated!
In fact, this rate is the same rate fully kinetic simulations which directly
solve the Vlasov equation obtain! "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "97232d57-a839-4b2e-8b22-4be31a7cb9ad",
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "# Compute the reconnection rate. Please note this may take several
minutes (5 - 7 minutes).\n",
+ "def compute_flux(pre, frame): \n",
+ " data_elc, data_ion, data_field = get_Gdata(pre, frame)\n",
+ " r_coords, r_fluid_elc, r_fluid_ion, r_field = read_data(data_elc,
data_ion, data_field)\n",
+ " # Fetch Ez at X point\n",
+ " EzX = r_field[int(Nx/2), int(Ny/2), 2]\n",
+ " # Compute vector potential at X point.\n",
+ " dx = r_coords[0][1] - r_coords[0][0]\n",
+ " dy = r_coords[1][1] - r_coords[1][0]\n",
+ " psi = np.zeros((Ny, Nx))\n",
+ " psi = calc_psi2d(r_field[...,3].transpose(),
r_field[...,4].transpose(), dx, dy)\n",
+ " psiX = psi[int(Ny/2), int(Nx/2)]\n",
+ " psi0 = np.max(np.abs(psi)) \n",
+ " return psiX, psi0, EzX\n",
+ " \n",
+ "num_frames = 61\n",
+ "psiX_5m = np.zeros(num_frames)\n",
+ "psi0_5m = np.zeros(num_frames)\n",
+ "EzX_5m = np.zeros(num_frames)\n",
+ "psiX_10m = np.zeros(num_frames)\n",
+ "psi0_10m = np.zeros(num_frames)\n",
+ "EzX_10m = np.zeros(num_frames)\n",
+ "t = np.linspace(0, 30/omega_ci, num_frames)\n",
+ "rate_fac = 1/(B0*vA_up)\n",
+ "for i in range(0, num_frames):\n",
+ " # Compute the flux at the X-point in resistive 5 moment simulation\n",
+ " psiX_5m[i], psi0_5m[i], EzX_5m[i] =
compute_flux(f\"{data_path}/resistive_gem_mass100/rt_5m_resistive_gem\", i) \n",
+ " # Compute the flux at the X-point in 10 moment simulation\n",
+ " psiX_10m[i], psi0_10m[i], EzX_10m[i] =
compute_flux(f\"{data_path}/10m_gem_mass100/rt_10m_gem\", i) \n",
+ "\n",
+ "# We can either use Ez or the time derivative of the vector potential psi
to estimate the reconnection rate. \n",
+ "psi_diff_5m = np.abs(psi0_5m - psiX_5m)\n",
+ "dpsi_dt_5m = np.gradient(psi_diff_5m, t)\n",
+ "psi_diff_10m = np.abs(psi0_10m - psiX_10m)\n",
+ "dpsi_dt_10m = np.gradient(psi_diff_10m, t)\n",
+ "plt.figure()\n",
+ "plt.plot(t*omega_ci, np.abs(dpsi_dt_5m)*rate_fac, 'k-', label=r'5m,
$d\\psi/dt$')\n",
+ "plt.plot(t*omega_ci, np.abs(EzX_5m)*rate_fac, 'k--', label=r'5m,
$E_z$')\n",
+ "plt.plot(t*omega_ci, np.abs(dpsi_dt_10m)*rate_fac, 'b-', label=r'10m,
$d\\psi/dt$')\n",
+ "plt.plot(t*omega_ci, np.abs(EzX_10m)*rate_fac, 'b--', label=r'10m,
$E_z$')\n",
+ "plt.xlabel(r'$t \\Omega_{ci}$')\n",
+ "plt.ylabel(r'$E_z/(v_A B_0)$')\n",
+ "plt.legend(loc='best')\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "02f2c461-511c-43a6-b24f-c4056e1fb481",
+ "metadata": {},
+ "source": [
+ "### Explaining \"fast\" reconnection with the Ohm's Law\n",
+ "\n",
+ "Let's now look at why the reconnection is faster in our multi-fluid
simulations by answering the following question: what is different about our
simulation compared to the Sweet-Parker argument? \n",
+ "1. MHD cannot get fast reconnection. We had to add at least two-fluid
effects to include dispersive waves. \n",
+ "2. Sweet-Parker assumed resistivity set the $\\delta/L$ scaling; that's
not true in our simulations. The layer width is bigger than if we only had
resistivity. \n",
+ "3. The layer is not genuinely incompressible. We need compressibility
so that the plasma can shock and push plasma in the outflow. \n",
+ "4. The conservation of energy argument is missing additional energy
flux terms which move energy into the layer.\n",
+ "\n",
+ "[Answer Form
5](https://docs.google.com/forms/d/e/1FAIpQLSc0le8Mn-4KTOqTZynZ9ztUfLtH7IhxW3mGi1x4SwURYmaWCQ/viewform?usp=dialog)\n",
+ "\n",
+ "Given this (key) difference between the Sweet-Parker argument and our
simulations, let's look at the effective Ohm's Law for our plasma. Unlike in
MHD, where the Ohm's Law prescribes the electric field,\n",
+ "\\begin{align}\n",
+ "\\mathbf{E}_{MHD} = \\mathbf{u} \\times \\mathbf{B} + \\eta
\\mathbf{J}\n",
+ "\\end{align}\n",
+ "we solved Maxwell's equations coupled to our multi-fluid system.
Nevertheless, we can still utilize the electron momentum equation to obtain an
Ohm's Law-like expression for the electric field:\n",
+ "\\begin{align}\n",
+ "\\mathbf{E} = \\frac{m_e}{q_e \\rho_e} \\left (\\frac{\\partial \\rho_e
\\mathbf{u}_e}{\\partial t} + \\nabla \\cdot \\left ( \\rho_e \\mathbf{u}_e
\\mathbf{u}_e + \\mathbf{P}_e \\right ) -\\frac{q_e}{m_e} \\rho_e \\mathbf{u}_e
\\times \\mathbf{B} - \\mathbf{R} \\right )\n",
+ "\\end{align}\n",
+ "which is often simplified using conservation of mass as \n",
+ "\\begin{align}\n",
+ "\\mathbf{E} = \\frac{m_e}{q_e} \\left (\\frac{\\partial
\\mathbf{u}_e}{\\partial t} + \\mathbf{u}_e \\cdot \\nabla \\mathbf{u}_e +
\\frac{1}{\\rho_e} \\nabla \\cdot \\mathbf{P}_e \\right ) -\\mathbf{u}_e
\\times \\mathbf{B} - \\frac{m_e}{q_e \\rho_e} \\mathbf{R}. \n",
+ "\\end{align}\n",
+ "\n",
+ "In the two-dimensional geometry of this magnetic reconnection setup, the
important electric field for supporting the reconnection is $E_z$. So our Ohm's
Law in general is\n",
+ "\\begin{align}\n",
+ "E_z = \\frac{m_e}{q_e} \\left (\\underbrace{\\frac{\\partial
u_{z_e}}{\\partial t} + u_{x_e} \\frac{\\partial u_{z_e}}{\\partial x} +
u_{y_e} \\frac{\\partial u_{z_e}}{\\partial y}}_{Inertia} + \\frac{1}{\\rho_e}
\\left [\\frac{\\partial P_{xz_e}}{\\partial x} + \\frac{\\partial
P_{yz_e}}{\\partial y} \\right ] \\right ) - \\underbrace{\\left (u_{x_e} B_y -
u_{y_e} B_x \\right )}_{Hall} - \\underbrace{\\frac{m_e}{q_e \\rho_e}
R_z}_{Resistivity}. \n",
+ "\\end{align}\n",
+ "with the off-diagonal pressure tensor components vanishing for the 5
moment simulations, and no resistive contribution to the 10 moment simulations.
\n",
+ "\n",
+ "Which of these terms can break field-lines? \n",
+ "1. Only resistivity. \n",
+ "2. Resistivity, inertia, and Hall term can break field lines (all the
terms in 5 moment). \n",
+ "3. Resistivity and the off-diagonal Pressure Tensor components. \n",
+ "4. All of the terms can break field lines.\n",
+ "\n",
+ "[Answer Form
6](https://docs.google.com/forms/d/e/1FAIpQLSfYwfxJJ8rMISttTfOlVlJEtzt4wIHQ4Z9QZHpHw_yA2FGR7A/viewform?usp=dialog)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "a085ca89-0379-4241-8c15-c375a567137e",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def fdx_2d_8th_order(data, dx, Nx, Ny):\n",
+ " diff_data = np.zeros((Nx, Ny))\n",
+ " for i in range(4, Nx-4):\n",
+ " for j in range(4, Ny-4): \n",
+ " diff_data[i,j] = (-data[i+4, j]/280.0 + data[i+3,
j]*4.0/105.0 - data[i+2, j]/5.0 + data[i+1, j]*4.0/5.0 \\\n",
+ " - data[i-1, j]*4.0/5.0 + data[i-2, j]/5.0 - data[i-3,
j]*4.0/105.0 + data[i-4, j]/280.0) / dx\n",
+ " return diff_data\n",
+ "\n",
+ "def fdy_2d_8th_order(data, dy, Nx, Ny):\n",
+ " diff_data = np.zeros((Nx, Ny))\n",
+ " for i in range(4, Nx-4):\n",
+ " for j in range(4, Ny-4): \n",
+ " diff_data[i,j] = (-data[i, j+4]/280.0 + data[i,
j+3]*4.0/105.0 - data[i, j+2]/5.0 + data[i, j+1]*4.0/5.0 \\\n",
+ " - data[i, j-1]*4.0/5.0 + data[i, j-2]/5.0 - data[i,
j-3]*4.0/105.0 + data[i, j-4]/280.0) / dy\n",
+ " return diff_data "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "61be64a7-bafa-4ed8-87de-93ffc390acde",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def ohm_10m(frame):\n",
+ " data_elc, data_ion, data_field =
get_Gdata(f\"{data_path}/10m_gem_mass100/rt_10m_gem\", frame)\n",
+ " r_coords, r_fluid_elc, r_fluid_ion, r_field = read_data(data_elc,
data_ion, data_field)\n",
+ " Nx = r_fluid_ion.shape[0]\n",
+ " Ny = r_fluid_ion.shape[1] \n",
+ " dx = r_coords[0][1] - r_coords[0][0]\n",
+ " dy = r_coords[1][1] - r_coords[1][0]\n",
+ " # Fetch mass density, in-plane momentum, out-of-plane velocity, and
pressure tensor Pij\n",
+ " rho_elc = r_fluid_elc[...,0]\n",
+ " rhoux_elc = r_fluid_elc[..., 1]\n",
+ " rhouy_elc = r_fluid_elc[..., 2]\n",
+ " _, uz_elc = pg.tools.get_vz(data_elc)\n",
+ " uz_elc = np.squeeze(uz_elc)\n",
+ " _, Pij_elc = pg.tools.get_pij(data_elc)\n",
+ " Pxz_elc = Pij_elc[...,2]\n",
+ " Pyz_elc = Pij_elc[...,4]\n",
+ " \n",
+ " # Compute 1/(q/m rho) div(P_iz) and 1/(q/m rho) (rhou_x du_z/dx +
rhou_y du_z/dy)\n",
+ " divRxz_elc = fdx_2d_8th_order(uz_elc, dx, Nx, Ny)\n",
+ " divRyz_elc = fdy_2d_8th_order(uz_elc, dy, Nx, Ny)\n",
+ " divPxz_elc = fdx_2d_8th_order(Pxz_elc, dx, Nx, Ny)\n",
+ " divPyz_elc = fdy_2d_8th_order(Pyz_elc, dy, Nx, Ny) \n",
+ "\n",
+ " divRxz_elc = rhoux_elc*divRxz_elc/(charge_elc/mass_elc*rho_elc)\n",
+ " divRyz_elc = rhouy_elc*divRyz_elc/(charge_elc/mass_elc*rho_elc) \n",
+ " divPxz_elc = divPxz_elc/(charge_elc/mass_elc*rho_elc)\n",
+ " divPyz_elc = divPyz_elc/(charge_elc/mass_elc*rho_elc) \n",
+ "\n",
+ " # Compute Hall term\n",
+ " Bx = r_field[...,3]\n",
+ " By = r_field[...,4]\n",
+ " hallz = np.zeros((Nx, Ny))\n",
+ " hallz = charge_elc/mass_elc*(rhoux_elc*By -
rhouy_elc*Bx)/(charge_elc/mass_elc*rho_elc)\n",
+ "\n",
+ " # Fetch Ez\n",
+ " Ez = r_field[...,2]\n",
+ "\n",
+ " return r_coords, Ez, hallz, divPxz_elc, divPyz_elc, divRxz_elc,
divRyz_elc"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d49ce79f-1e2e-45ce-b58a-d640f6700a45",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def ohm_5m(frame):\n",
+ " data_elc, data_ion, data_field =
get_Gdata(f\"{data_path}/resistive_gem_mass100/rt_5m_resistive_gem\", frame)\n",
+ " r_coords, r_fluid_elc, r_fluid_ion, r_field = read_data(data_elc,
data_ion, data_field)\n",
+ " Nx = r_fluid_ion.shape[0]\n",
+ " Ny = r_fluid_ion.shape[1] \n",
+ " dx = r_coords[0][1] - r_coords[0][0]\n",
+ " dy = r_coords[1][1] - r_coords[1][0]\n",
+ " # Compute Hall term\n",
+ " rho_elc = r_fluid_elc[...,0]\n",
+ " rho_ion = r_fluid_ion[...,0]\n",
+ " rhoux_elc = r_fluid_elc[..., 1]\n",
+ " rhouy_elc = r_fluid_elc[..., 2]\n",
+ " Bx = r_field[...,3]\n",
+ " By = r_field[...,4]\n",
+ " hallz = np.zeros((Nx, Ny))\n",
+ " hallz = charge_elc/mass_elc*(rhoux_elc*By -
rhouy_elc*Bx)/(charge_elc/mass_elc*rho_elc)\n",
+ "\n",
+ " # Fetch Ez\n",
+ " Ez = r_field[...,2]\n",
+ "\n",
+ " # Compute resistivity\n",
+ " _, uz_elc = pg.tools.get_vz(data_elc)\n",
+ " _, uz_ion = pg.tools.get_vz(data_ion)\n",
+ " _, T_elc = pg.tools.get_temp(data_elc, gas_gamma)\n",
+ " uz_elc = np.squeeze(uz_elc)\n",
+ " uz_ion = np.squeeze(uz_ion)\n",
+ " T_elc_norm = np.squeeze(T_elc/T_elc_ref)\n",
+ " tau_ei =
Lambda_ee/np.log(Lambda_ee)*6.0*np.sqrt(2.0*np.pi**3*T_elc_norm**3*mass_elc*mass_ion**2)/rho_ion\n",
+ " eta_z = (-rho_elc/tau_ei*(uz_elc -
uz_ion))/(charge_elc/mass_elc*rho_elc)\n",
+ "\n",
+ " # Compute 1/(q/m rho) (rhou_x du_z/dx + rhou_y du_z/dy)\n",
+ " divRxz_elc = fdx_2d_8th_order(uz_elc, dx, Nx, Ny)\n",
+ " divRyz_elc = fdy_2d_8th_order(uz_elc, dy, Nx, Ny)\n",
+ "\n",
+ " divRxz_elc = rhoux_elc*divRxz_elc/(charge_elc/mass_elc*rho_elc)\n",
+ " divRyz_elc = rhouy_elc*divRyz_elc/(charge_elc/mass_elc*rho_elc) \n",
+ " \n",
+ " return r_coords, Ez, hallz, eta_z, divRxz_elc, divRyz_elc"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "ef3518a8-8b2a-49d0-8ed1-377ac2430e70",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Choose time = 20 Omega_{ci}^{-1} near the peak of the reconnection
(frame 40)\n",
+ "coords, Ez, hallz, eta_z, divRxz_elc, divRyz_elc = ohm_5m(40)\n",
+ "sum_ohm = -hallz - eta_z + divRxz_elc + divRyz_elc\n",
+ "\n",
+ "x_point = int(Nx/2)\n",
+ "plt.figure()\n",
+ "plt.plot(coords[1]/di, Ez[x_point, :], 'k', label=r'$E_z$')\n",
+ "plt.plot(coords[1]/di, -hallz[x_point, :], 'r', label=r'$-u_e x B$')\n",
+ "plt.plot(coords[1]/di, -eta_z[x_point, :], 'b', label=r'$-1/(qn) \\eta
J_z$')\n",
+ "plt.plot(coords[1]/di, divRxz_elc[x_point, :], 'y', label=r'$u_x
\\partial_x u_z$')\n",
+ "plt.plot(coords[1]/di, divRyz_elc[x_point, :], 'c', label=r'$u_y
\\partial_y u_z$')\n",
+ "plt.plot(coords[1]/di, sum_ohm[x_point, :], 'm', label=r'Sum')\n",
+ "plt.xlabel(r'$y/d_i$')\n",
+ "plt.legend(loc='upper right')\n",
+ "plt.xlim(-6.0, 6.0)\n",
+ "\n",
+ "plt.figure()\n",
+ "plt.plot(coords[1]/di, Ez[x_point, :], 'k', label=r'$E_z$')\n",
+ "plt.plot(coords[1]/di, -hallz[x_point, :], 'r', label=r'$-u_e x B$')\n",
+ "plt.plot(coords[1]/di, -eta_z[x_point, :], 'b', label=r'$-1/(qn) \\eta
J_z$')\n",
+ "plt.plot(coords[1]/di, divRxz_elc[x_point, :], 'y', label=r'$u_x
\\partial_x u_z$')\n",
+ "plt.plot(coords[1]/di, divRyz_elc[x_point, :], 'c', label=r'$u_y
\\partial_y u_z$')\n",
+ "plt.plot(coords[1]/di, sum_ohm[x_point, :], 'm', label=r'Sum')\n",
+ "plt.xlabel(r'$y/d_i$')\n",
+ "plt.legend(loc='upper right')\n",
+ "plt.xlim(-1.0, 1.0)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "465609cc-8e3e-4f84-af19-08f5acc65e27",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Choose time = 20 Omega_{ci}^{-1} near the peak of the reconnection
(frame 40)\n",
+ "coords, Ez, hallz, divPxz_elc, divPyz_elc, divRxz_elc, divRyz_elc =
ohm_10m(40)\n",
+ "sum_ohm = -hallz + divPxz_elc + divPyz_elc + divRxz_elc + divRyz_elc\n",
+ "\n",
+ "x_point = int(Nx/2)\n",
+ "plt.figure()\n",
+ "plt.plot(coords[1]/di, Ez[x_point, :], 'k', label=r'$E_z$')\n",
+ "plt.plot(coords[1]/di, -hallz[x_point, :], 'r', label=r'$-u_e x B$')\n",
+ "plt.plot(coords[1]/di, divPxz_elc[x_point, :], 'b', label=r'$1/(qn)
\\partial_x P_{xz}$')\n",
+ "plt.plot(coords[1]/di, divPyz_elc[x_point, :], 'g', label=r'$1/(qn)
\\partial_y P_{yz}$')\n",
+ "plt.plot(coords[1]/di, divRxz_elc[x_point, :], 'y', label=r'$u_x
\\partial_y u_z$')\n",
+ "plt.plot(coords[1]/di, divRyz_elc[x_point, :], 'c', label=r'$u_y
\\partial_y u_z$')\n",
+ "plt.plot(coords[1]/di, sum_ohm[x_point, :], 'm', label=r'Sum')\n",
+ "plt.xlabel(r'$y/d_i$')\n",
+ "plt.legend(loc='upper right')\n",
+ "plt.xlim(-6.0, 6.0)\n",
+ "\n",
+ "plt.figure()\n",
+ "plt.plot(coords[1]/di, Ez[x_point, :], 'k', label=r'$E_z$')\n",
+ "plt.plot(coords[1]/di, -hallz[x_point, :], 'r', label=r'$-u_e x B$')\n",
+ "plt.plot(coords[1]/di, divPxz_elc[x_point, :], 'b', label=r'$1/(qn)
\\partial_x P_{xz}$')\n",
+ "plt.plot(coords[1]/di, divPyz_elc[x_point, :], 'g', label=r'$1/(qn)
\\partial_y P_{yz}$')\n",
+ "plt.plot(coords[1]/di, divRxz_elc[x_point, :], 'y', label=r'$u_x
\\partial_x u_z$')\n",
+ "plt.plot(coords[1]/di, divRyz_elc[x_point, :], 'c', label=r'$u_y
\\partial_y u_z$')\n",
+ "plt.plot(coords[1]/di, sum_ohm[x_point, :], 'm', label=r'Sum')\n",
+ "plt.xlabel(r'$y/d_i$')\n",
+ "plt.legend(loc='upper right')\n",
+ "plt.xlim(-1.0, 1.0)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "602b46ec-1aac-4e8c-8034-14381daeef13",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "%stop_runtime test_cpu\n",
+ "%switch_runtime local"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1c77260d-6ed4-4cd1-be96-099fe826627a",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.11.6"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git
a/dev-tools/deployment/jupyterhub/user-container/MD/data/Gkeyll/plotE_z.ipynb
b/dev-tools/deployment/jupyterhub/user-container/MD/data/Gkeyll/post_gkeyll_plotE_z.ipynb
similarity index 71%
rename from
dev-tools/deployment/jupyterhub/user-container/MD/data/Gkeyll/plotE_z.ipynb
rename to
dev-tools/deployment/jupyterhub/user-container/MD/data/Gkeyll/post_gkeyll_plotE_z.ipynb
index ff7cc617b0..e09e0fe5c7 100644
---
a/dev-tools/deployment/jupyterhub/user-container/MD/data/Gkeyll/plotE_z.ipynb
+++
b/dev-tools/deployment/jupyterhub/user-container/MD/data/Gkeyll/post_gkeyll_plotE_z.ipynb
@@ -2,15 +2,19 @@
"cells": [
{
"cell_type": "markdown",
- "id": "ec13f649-d9ee-47a2-9761-131abb1e2858",
+ "id": "d44ab5e6-1229-438e-a780-a73cf2855abf",
"metadata": {},
"source": [
"#### Initialize, Authenticate, and Configure Remote Execution
Environment\n",
"\n",
"##### Start a Cybershuttle Remote Executor: This will submit a remote job
to Anvil\n",
- "- **Cluster** - Anvil\n",
+ "- **Cluster** - Anvil & Jetstream\n",
"- **Community Allocation** - Gkeyll\n",
- "- **Requested Resources** - 4 CPUs, 4GB Memory, 60 Minutes, Shared Queue"
+ "- **Requested Resources** \n",
+ " - 4 CPUs, 4GB Memory, 60 Minutes, Shared Queue\n",
+ " - 1 CPUs, 60 Minutes, Cloud Queue\n",
+ "\n",
+ "##### This remote job enables the generation of plotE_z with post gkeyll
code.\n"
]
},
{
@@ -20,12 +24,21 @@
"metadata": {},
"outputs": [],
"source": [
- "%pip install --force-reinstall airavata-jupyter-magic\n",
+ "%pip install --force-reinstall airavata-jupyter-magic > /dev/null\n",
"import airavata_jupyter_magic\n",
"\n",
"%authenticate\n",
+ "\n",
+ "# PLEASE NOTE: At a given time, ONLY run a single job in one of the
available clusters: Anvil OR Jetstream. \n",
+ "#Anvil is the default cluster and if you need to run on Jetstream;
comment the Anvil and uncomment the Jetstream\n",
+ "\n",
+ "#Running on Anvil remote cluster\n",
"%request_runtime test_cpu --cluster=Anvil --cpus=4 --memory=4096
--walltime=60 --queue=shared --group=Gkeyll\n",
- "%switch_runtime test_cpu"
+ "%switch_runtime test_cpu\n",
+ "\n",
+ "#Running on Jetstream cluster\n",
+ "#%request_runtime test_cpu --cluster=JS --cpus=1 --walltime=60
--queue=cloud --group=Gkeyll\n",
+ "#%switch_runtime test_cpu"
]
},
{
@@ -35,7 +48,15 @@
"metadata": {},
"outputs": [],
"source": [
- "!ls /anvil/projects/x-phy220105/gkylMarch2025/vlabApps/gemReconnection"
+ "#Defininf the data paths in remote clusters\n",
+ "%ls -l /data\n",
+ "\n",
+ "#Anvil and Jetstream data path\n",
+ "data_path = \"/data/gemReconnection\"\n",
+ "\n",
+ "!ls /data/gemReconnection\n",
+ "\n",
+ "#Note: This cell may execute longer if the job is not yet active in the
remote cluster. Please pause till the files are listed below"
]
},
{
@@ -49,14 +70,13 @@
"import matplotlib.pyplot as plt\n",
"import postgkyl as pg\n",
"\n",
- "data_base_path =
\"/anvil/projects/x-phy220105/gkylMarch2025/vlabApps/gemReconnection\"\n",
"run = Path.cwd()\n",
"\n",
"def getModelType():\n",
" frame = 0\n",
" models = [\"5m\", \"10m\"]\n",
" for model in models:\n",
- " path =
Path(f\"{data_base_path}/rt_{model}_gem_gzero-field_{frame}.gkyl\")\n",
+ " path =
Path(f\"{data_path}/rt_{model}_gem_gzero-field_{frame}.gkyl\")\n",
" if path.is_file():\n",
" return model\n",
" error = \"Failed to find input file \" + str(path)\n",
@@ -72,7 +92,7 @@
"source": [
"frame = 0\n",
"model = getModelType()\n",
- "filename = run /
f\"{data_base_path}/rt_{model}_gem_gzero-field_{frame}.gkyl\"\n",
+ "filename = run /
f\"{data_path}/rt_{model}_gem_gzero-field_{frame}.gkyl\"\n",
"filename = str(filename)\n",
"print(filename)"
]
diff --git a/dev-tools/deployment/jupyterhub/user-container/MD/init.sh
b/dev-tools/deployment/jupyterhub/user-container/MD/init.sh
index 2494ef4249..2882007d53 100755
--- a/dev-tools/deployment/jupyterhub/user-container/MD/init.sh
+++ b/dev-tools/deployment/jupyterhub/user-container/MD/init.sh
@@ -1,8 +1,6 @@
#!/bin/bash
TARGET_DIR="/home/jovyan/work"
-SHARED_TMP="/home/jovyan/shared_data_tmp"
-
mkdir -p "$TARGET_DIR"
@@ -16,12 +14,4 @@ else
echo "Docker default files already exist, skipping copy."
fi
-
-if [ -d "$SHARED_TMP" ]; then
- echo "Creating symlinks for shared read-only data inside $TARGET_DIR..."
- for item in "$SHARED_TMP"/*; do
- ln -s "$item" "$TARGET_DIR/$(basename "$item")"
- done
-fi
-
exec "$@"
