Prepared uniform plasma case for Cherenkov test

3 files changed, 218 insertions, 23 deletions

diff --git a/Exec/uniform_plasma/ParticleProb.cpp b/Exec/uniform_plasma/ParticleProb.cpp
index e23be8f0b..cef618f61 100644
--- a/Exec/uniform_plasma/ParticleProb.cpp
+++ b/Exec/uniform_plasma/ParticleProb.cpp
@@ -19,15 +19,23 @@ PhysicalParticleContainer::InitData()
 {
     BL_PROFILE("PPC::InitData()");
 
-    charge = -PhysConst::q_e;
-    mass = PhysConst::m_e;
+    // species_id 0 : electrons
+    // species_id 1 : Hydrogen ions
+    // Note: the ions of the plasma are implicitly motionless, and so are not part of the simulation
+    if (species_id == 0) {
+        charge = -PhysConst::q_e;
+    	mass = PhysConst::m_e;
+    } else {
+        charge = PhysConst::q_e;
+        mass = PhysConst::m_p;
+    }
 
     const int lev = 0;
 
     const Geometry& geom = Geom(lev);
     const Real* dx  = geom.CellSize();
 
-    Real weight, u_th;
+    Real weight, u_th, ux_m, uy_m, uz_m;
     Real particle_shift_x, particle_shift_y, particle_shift_z;
     int n_part_per_cell;
     {
@@ -43,8 +51,17 @@ PhysicalParticleContainer::InitData()
       #endif
 
       u_th = 0.;
+      ux_m = 0.;
+      uy_m = 0.;
+      uz_m = 0.;
       pp.query("u_th", u_th);
+      pp.query("ux_m", ux_m);
+      pp.query("uy_m", uy_m);
+      pp.query("uz_m", uz_m);
       u_th *= PhysConst::c;
+      ux_m *= PhysConst::c;
+      uy_m *= PhysConst::c;
+      uz_m *= PhysConst::c;
     }
 
     std::array<Real,PIdx::nattribs> attribs;
@@ -54,7 +71,6 @@ PhysicalParticleContainer::InitData()
     // Initialize random generator for normal distribution
     std::default_random_engine generator;
     std::normal_distribution<Real> velocity_distribution(0.0, u_th);
-    std::uniform_real_distribution<Real> position_distribution(0.0,1.0);
 
     for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi)
     {
@@ -70,28 +86,26 @@ PhysicalParticleContainer::InitData()
             for (int i_part=0; i_part<n_part_per_cell;i_part++)
             {
                 // Randomly generate the speed according to a normal distribution
-                Real ux = velocity_distribution(generator);
-                Real uy = velocity_distribution(generator);
-                Real uz = velocity_distribution(generator);
+                Real ux_th = velocity_distribution(generator);
+                Real uy_th = velocity_distribution(generator);
+                Real uz_th = velocity_distribution(generator);
 
                 // Randomly generate the positions (uniformly inside each cell)
-                Real particle_shift_x = position_distribution(generator);
-                Real particle_shift_y = position_distribution(generator);
-                Real particle_shift_z = position_distribution(generator);
+                Real particle_shift = (0.5+i_part)/n_part_per_cell;
 
 #if (BL_SPACEDIM == 3)
-                Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift_x)*dx[0];
-                Real y = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift_y)*dx[1];
-                Real z = tile_real_box.lo(2) + (iv[2]-boxlo[2] + particle_shift_z)*dx[2];
+                Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0];
+                Real y = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1];
+                Real z = tile_real_box.lo(2) + (iv[2]-boxlo[2] + particle_shift)*dx[2];
 #elif (BL_SPACEDIM == 2)
-                Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift_x)*dx[0];
+                Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0];
                 Real y = 0.0;
-                Real z = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift_z)*dx[1];
-#endif   
+                Real z = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1];
+#endif
 
-                attribs[PIdx::ux] = ux;
-                attribs[PIdx::uy] = uy;
-                attribs[PIdx::uz] = uz;
+                attribs[PIdx::ux] = ux_m + ux_th;
+                attribs[PIdx::uy] = uy_m + uy_th;
+                attribs[PIdx::uz] = uz_m + uz_th;
 
                 AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs);
             }
diff --git a/Exec/uniform_plasma/Visualization.ipynb b/Exec/uniform_plasma/Visualization.ipynb
new file mode 100644
index 000000000..559f36482
--- /dev/null
+++ b/Exec/uniform_plasma/Visualization.ipynb
@@ -0,0 +1,177 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "deletable": true,
+    "editable": true
+   },
+   "source": [
+    "# Overview\n",
+    "\n",
+    "This a notebook that inspects the results of a WarpX simulation.\n",
+    "\n",
+    "# Instructions\n",
+    "\n",
+    "Execute the cells below one by one, by selecting them with your mouse and typing `Shift + Enter`"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "deletable": true,
+    "editable": true
+   },
+   "outputs": [],
+   "source": [
+    "# General imports\n",
+    "import glob, sys\n",
+    "# Import numpy and matplotlib\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "%matplotlib\n",
+    "# Import yt\n",
+    "import yt\n",
+    "yt.funcs.mylog.setLevel(50) # Deactivate yt verbose\n",
+    "\n",
+    "# Import interactive notebook feature\n",
+    "#from IPython.display import clear_output\n",
+    "from ipywidgets import interact, IntSlider, RadioButtons\n",
+    "\n",
+    "# Find output iterations\n",
+    "file_list = glob.glob('plt?????')\n",
+    "iterations = [ int(file_name[3:]) for file_name in file_list ]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "deletable": true,
+    "editable": true
+   },
+   "source": [
+    "# Functions to extract and plot the fields"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "def extra_data( field, iteration ):\n",
+    "    \"\"\"\n",
+    "    Extract the 3D array that corresponds to `field` at the \n",
+    "    requested iteration, and returns as a `YTarray`, which \n",
+    "    can be used exactly like a numpy array.\n",
+    "    \"\"\"\n",
+    "    ds = yt.load( './plt%05d/' %iteration )\n",
+    "    data = ds.covering_grid(level=0, \n",
+    "            left_edge=ds.domain_left_edge, \n",
+    "            dims=ds.domain_dimensions)\n",
+    "    return( data[field] ) "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "deletable": true,
+    "editable": true
+   },
+   "outputs": [],
+   "source": [
+    "def plot_field( iteration, field, slicing_direction='y' ):\n",
+    "    # Get the data\n",
+    "    field_array = extra_data( field, iteration )\n",
+    "    \n",
+    "    # Get the position of the edges\n",
+    "    ds = yt.load( './plt%05d/' %iteration )\n",
+    "    left_edge = ds.domain_left_edge.convert_to_mks()*1.e6\n",
+    "    right_edge = ds.domain_right_edge.convert_to_mks()*1.e6\n",
+    "\n",
+    "    if slicing_direction == 'x':\n",
+    "        n = int( ds.domain_dimensions[0]//2 )\n",
+    "        data2d = field_array[n, :, :]\n",
+    "        extent = [ left_edge[2], right_edge[2], left_edge[1], right_edge[1] ]\n",
+    "    elif slicing_direction == 'y':\n",
+    "        n = int( ds.domain_dimensions[1]//2 )\n",
+    "        data2d = field_array[:, n, :]\n",
+    "        extent = [ left_edge[2], right_edge[2], left_edge[0], right_edge[0] ]\n",
+    "    elif slicing_direction == 'z':\n",
+    "        n = int( ds.domain_dimensions[2]//2 )\n",
+    "        data2d = field_array[:, :, n]\n",
+    "        extent = [ left_edge[1], right_edge[1], left_edge[0], right_edge[0] ]\n",
+    "    plt.clf()\n",
+    "    plt.imshow( data2d, interpolation='nearest', cmap='viridis',\n",
+    "              origin='lower', extent=extent )\n",
+    "    plt.colorbar()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "deletable": true,
+    "editable": true
+   },
+   "source": [
+    "# Interactive viewer"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "deletable": true,
+    "editable": true
+   },
+   "outputs": [],
+   "source": [
+    "interact(plot_field, \n",
+    "         iteration = IntSlider(min=min(iterations), max=max(iterations), step=iterations[1]-iterations[0]),\n",
+    "         field = RadioButtons( options=['jx', 'jy', 'jz', 'Ex', 'Ey', 'Ez'], value='jz'),\n",
+    "         slicing_direction = RadioButtons( options=[ 'x', 'y', 'z'], value='y') )"
+   ]
+  }
+ ],
+ "metadata": {
+  "anaconda-cloud": {},
+  "kernelspec": {
+   "display_name": "Python 3",
+   "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.6.0"
+  },
+  "widgets": {
+   "state": {
+    "5ca3ffc563434f7487da7b1bd657761b": {
+     "views": [
+      {
+       "cell_index": 5
+      }
+     ]
+    }
+   },
+   "version": "1.2.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 1
+}
diff --git a/Exec/uniform_plasma/inputs b/Exec/uniform_plasma/inputs
index 21eb0c7db..5e87c689e 100644
--- a/Exec/uniform_plasma/inputs
+++ b/Exec/uniform_plasma/inputs
@@ -1,19 +1,22 @@
 # Maximum number of time steps
-max_step = 20
+max_step = 100
 
 # number of grid points
-amr.n_cell =  128 128 128
+amr.n_cell =  64 64 64
 
-# Maximum allowable size of each subdomain in the problem domain; 
+amr.plot_int = 5   # How often to write plotfiles.
+
+# Maximum allowable size of each subdomain in the problem domain;
 #    this is used to decompose the domain for parallel calculations.
 amr.max_grid_size = 32
 
 # Maximum level in hierarchy (for now must be 0, i.e., one level in total)
 amr.max_level = 0
+particles.nspecies = 2
 
 # Geometry
 geometry.coord_sys   = 0                  # 0: Cartesian
-geometry.is_periodic = 1     1     1      # Is periodic?  
+geometry.is_periodic = 1     1     1      # Is periodic?
 geometry.prob_lo     = -20.e-6   -20.e-6   -20.e-6    # physical domain
 geometry.prob_hi     =  20.e-6    20.e-6    20.e-6
 
@@ -32,3 +35,4 @@ warpx.cfl = 1.0
 uniform_plasma.num_particles_per_cell = 40
 uniform_plasma.n_e = 1.e25  # number of electrons per m^3
 uniform_plasma.u_th  = 0.01 # uth the std of the (unitless) momentum
+uniform_plasma.uz_m  = 10.  # Mean momentum along x (unitless)