Testing plots and new collision operator in documentation

docs/bibs/fokkerplanck.bib

@@ -11,3 +11,16 @@ @article{Lenard1958
 year = {1958}
 }
 
+@article{Dougherty1964,
+abstract = {A model Fokker-Planck equation for a single species of particles in a plasma is discussed. This equation has several properties in common with the real equation, and ascribes an approximately correct behavior to most of the particles, though incorrect for the high-energy tail in the thermal distribution. It is shown that the equation can be solved completely, for the small perturbations of a uniform plasma by electric fields harmonic in space and time, in an external magnetic field. Applications to ionospheric radar scattering are briefly discussed; it is shown that in certain circumstances, ion-ion collisions can have a profound effect on the scattering even though the collision frequency is much smaller than the ion gyrofrequency, and this appears to agree with observation.},
+author = {Dougherty, J. P.},
+doi = {10.1063/1.2746779},
+issn = {00319171},
+journal = {Physics of Fluids},
+number = {11},
+pages = {1788},
+title = {{Model Fokker-Planck Equation for a Plasma and Its Solution}},
+url = {https://aip.scitation.org/doi/10.1063/1.2746779},
+volume = {7},
+year = {1964}
+}

docs/bibs/usage_details.bib

@@ -0,0 +1,15 @@
+@article{Joglekar2018,
+author = {Joglekar, Archis S and Winjum, Benjamin J and Tableman, Adam and Wen, Han and Tzoufras, Michail and Mori, Warren B},
+doi = {10.1088/1361-6587/aab978},
+issn = {0741-3335},
+journal = {Plasma Physics and Controlled Fusion},
+keywords = {and length-scales as well,appear in colour only,as plasma densities and,cartesian tensor,fokker,in the online journal,kinetic modeling,occur over disparate time-,planck,plasma physics phenomenon can,some fi gures may,spherical harmonic,tempera-,vlasov},
+month = {jun},
+number = {6},
+pages = {064010},
+publisher = {IOP Publishing},
+title = {{Validation of OSHUN against collisionless and collisional plasma physics}},
+url = {http://iopscience.iop.org/article/10.1088/1361-6587/aab978 http://stacks.iop.org/0741-3335/60/i=6/a=064010?key=crossref.71a4d809c09e954863819a819c5ceb29},
+volume = {60},
+year = {2018}
+}

docs/electricfield.rst

@@ -20,4 +20,8 @@ given by
 
 This solver is tested to reproduce analytical solutions to a periodic Poisson system.
 
+The Poisson Solver is tested in `tests/test_fieldsolver.py` and the tests are illustrated in
+`notebooks/test_poisson.ipynb`. Below, we show a sample test result.
 
+.. image:: images/poisson_solver.png
+   :width: 900

docs/fokker-planck.rst

@@ -1,27 +1,32 @@
 Solving the Fokker-Planck Equation
 ----------------------------------------
 
-There are many approximations to the Fokker-Planck equation as applied to a distribution function of electrons.
-
-Lenard-Bernstein Form
-****************************
-VlaPy contains the Lenard-Bernstein :cite:`Lenard1958` implementation of the Fokker-Planck equation. This equation is given by
+VlaPy implements the Lenard-Bernstein :cite:`Lenard1958` (LB) approximation and the Dougherty :cite:`Dougherty1964` (DG)
+approximation of the Fokker-Planck equation. The LB approximation is given by
 
 .. math::
     \frac{\partial f}{\partial t} = \nu_{ee} \frac{\partial}{\partial v} \left(v f + v_0^2 \frac{\partial f}{\partial v} \right)
 
 where :math:`\nu_{ee}` and :math:`v_0` are the electron-electron collision frequency, and thermal velocity, respectively.
 
+The DG approximation is given by
+
+.. math::
+    \frac{\partial f}{\partial t} = \nu_{ee} \frac{\partial}{\partial v} \left((v - \underline{v}) f + v_t^2 \frac{\partial f}{\partial v} \right)
+
+where :math:`\underline{v} = \int f v dv` and :math:`v_t^2 = \int f (v - \underline{v})^2 dv` are the mean electron
+velocity, and the shifted thermal electron velocity.
+
 
 Differencing Scheme
 ====================
 
 This operator is differenced backwards in time, and center differenced in velocity space, which gives
 
 .. math::
-    \frac{f^{n+1}_{\alpha} - f^{n}_{\alpha}}{\Delta t} = \nu_{ee} \left[f^{n+1}_\alpha + v_\alpha \Delta_v(f^{n+1}_{\alpha}) + v_0^2 \Delta^2_v(f^{n+1}_{\alpha})\right]
+    \frac{f^{n+1}_{\alpha} - f^{n}_{\alpha}}{\Delta t} = \nu_{ee} \left[f^{n+1}_\alpha + \bar{v}_\alpha \Delta_v(f^{n+1}_{\alpha}) + v_{rms}^2 \Delta^2_v(f^{n+1}_{\alpha})\right]
 
-where
+where :math:`\bar{v} = v, v_{rms}^2 = \int f v^2 dv` for the LB operator and :math:`\bar{v} = v - \underline{v}, v_{rms}^2 = \int f \bar{v}^2 dv` for the DG operator.
 
 .. math::
     \Delta_v(f^{n+1}_{\alpha})= \frac{f^{n+1}_{\alpha+1} - f^{n+1}_{\alpha-1}}{2\Delta v}
@@ -46,9 +51,22 @@ Tests
 
 This solver is tested to
 1) return df/dt = 0 if a Maxwell-Boltzmann distribution is provided as input
-2) conserve energy and density
-3) relax to a Maxwellian of the right temperature and without a drift velocity
+2) conserve density, energy, and depending on the operator, velocity.
+3) relax to an operator-dependent Maxwellian of the right temperature and drift velocity.
+
+These tests are illustrated in `notebooks/test_fokker_planck.ipynb` and below:
+
+.. image:: images/Maxwell_Solution.png
+   :width: 900
+
+.. image:: images/LB_conservation.png
+   :width: 900
+
+.. image:: images/LB_no_conservation.png
+   :width: 900
 
+.. image:: images/DG_conservation.png
+   :width: 900
 
 
 .. bibliography:: bibs/fokkerplanck.bib

docs/usage_details.rst

@@ -15,36 +15,54 @@ For example, `run_vlapy.py` has the following code
 
 .. code-block:: python
 
+    all_params_dict = {
+        "nx": 48,
+        "xmin": 0.0,
+        "xmax": 2.0 * np.pi / 0.3,
+        "nv": 512,
+        "vmax": 6.0,
+        "nt": 1000,
+        "tmax": 100,
+        "nu": 0.0,
+    }
+
+    pulse_dictionary = {
+        "first pulse": {
+            "start_time": 0,
+            "rise_time": 5,
+            "flat_time": 10,
+            "fall_time": 5,
+            "a0": 1e-6,
+            "k0": 0.3,
+        }
+    }
+
+    params_to_log = ["w0", "k0", "a0"]
+
+    pulse_dictionary["first pulse"]["w0"] = np.real(
+        z_function.get_roots_to_electrostatic_dispersion(
+            wp_e=1.0, vth_e=1.0, k0=pulse_dictionary["first pulse"]["k0"]
+        )
+    )
+
+    mlflow_exp_name = "Landau Damping-test"
+
     manager.start_run(
-        nx=48,
-        nv=512,
-        nt=1000,
-        tmax=100,
-        nu=0.001,
-        w0=1.1598,
-        k0=0.3,
-        a0=1e-5,
-        diagnostics=landau_damping.LandauDamping(),
-        name="Landau Damping",
+        all_params=all_params_dict,
+        pulse_dictionary=pulse_dictionary,
+        diagnostics=landau_damping.LandauDamping(params_to_log),
+        name=mlflow_exp_name,
     )
 
-Through this interface, the user can specify the boundaries of the domain in space and time. The user can
-also specify initial conditions such as the settings of the wave-driver module, the collision frequency `nu`, and custom
+Through this interface, the user can specify the initial conditions of the domain in space and time. The user can
+also specify the settings of the wave-driver module, the collision frequency `nu`, and custom
 diagnostic routines.
 
 
 Wave Driver
 ===============
-The wave driver module currently has a fixed form given by
-
-.. code-block:: python
-
-    def driver_function(x, t):
-        envelope = np.exp(-((t - 8) ** 8.0) / 4.0 ** 8.0)
-        return envelope * a0 * np.cos(k0 * x + w0 * t)
-
-This effectively mimics a flat-top pulse in time. The ``w0, a0,`` and ``k0`` quantities can be specified to generate waves
-of specified frequency, amplitude, and wave-number, respectively.
+The wave driver module uses a 5th order polynomial that is given in ref. [@Joglekar2018]. The implementation
+can be found in `vlapy/manager.py`.
 
 
 Simulation Management

paper/paper.md

@@ -163,4 +163,6 @@ We use xarray [@Hoyer2017] for file storage and MLFlow [@Zaharia2018] for experi
 
 We acknowledge valuable discussions with Pierre Navarro on the implementation of the Vlasov equation.
 
+We are grateful for the editors' and reviewers' thorough feedback that improved the software as well as manuscript.
+
 # References