From 5033acd8f366f9917945847acb03087035cfa685 Mon Sep 17 00:00:00 2001
From: Casey Icenhour <cticenho@ncsu.edu>
Date: Thu, 13 Sep 2018 11:15:30 -0600
Subject: [PATCH] Fixes for EnergyFluxBC (residual), EnergyExcitation and
 Energy Ionization (Jacobians), and JouleHeating (OffDiag Jacobians)

---
 src/bcs/EnergyFluxBC.C         | 9 ++++-----
 src/kernels/EnergyExcitation.C | 6 +++++-
 src/kernels/EnergyIonization.C | 6 +++++-
 src/kernels/JouleHeating.C     | 5 ++---
 4 files changed, 16 insertions(+), 10 deletions(-)

diff --git a/src/bcs/EnergyFluxBC.C b/src/bcs/EnergyFluxBC.C
index 1b55e05..eebbf7e 100644
--- a/src/bcs/EnergyFluxBC.C
+++ b/src/bcs/EnergyFluxBC.C
@@ -40,7 +40,7 @@ EnergyFluxBC::computeQpResidual()
 {
   return _test[_i][_qp] * _normals[_qp] *
          ((5 / 2) * (2 / 3) * _u[_qp] * 0.25 *
-              std::pow((8 * (2 / 3) * _u[_qp]) / (libMesh::pi * _electron_mass), 0.25) *
+              std::pow((8 * (2 / 3) * _u[_qp]) / (libMesh::pi * _electron_mass), 0.5) *
               _electron_density[_qp] * _normals[_qp] -
           _secondary_electron_coefficient * (5 / 2) * _T_secondary_electron *
               (-_ion_mobility * _ion_density[_qp] * _grad_potential[_qp]));
@@ -50,9 +50,8 @@ Real
 EnergyFluxBC::computeQpJacobian()
 {
   return _test[_i][_qp] * _normals[_qp] *
-         ((5 / 2) * (2 / 3) * 0.25 *
-          std::pow((8 * (2 / 3)) / (libMesh::pi * _electron_mass), 0.25) * _electron_density[_qp] *
-          1.5 * std::pow(_u[_qp], 0.25) * _phi[_j][_qp] * _normals[_qp]);
+         ((5 / 2) * (2 / 3) * 0.25 * std::pow((8 * (2 / 3)) / (libMesh::pi * _electron_mass), 0.5) *
+          _electron_density[_qp] * 1.5 * std::pow(_u[_qp], 0.5) * _phi[_j][_qp] * _normals[_qp]);
 }
 
 Real
@@ -61,7 +60,7 @@ EnergyFluxBC::computeQpOffDiagJacobian(unsigned int jvar)
   if (jvar == _electron_id)
     return _test[_i][_qp] * _normals[_qp] *
            ((5 / 2) * (2 / 3) * _u[_qp] * 0.25 *
-            std::pow((8 * (2 / 3) * _u[_qp]) / (libMesh::pi * _electron_mass), 0.25) *
+            std::pow((8 * (2 / 3) * _u[_qp]) / (libMesh::pi * _electron_mass), 0.5) *
             _phi[_j][_qp] * _normals[_qp]);
 
   else if (jvar == _ion_id)
diff --git a/src/kernels/EnergyExcitation.C b/src/kernels/EnergyExcitation.C
index 82d4e19..4dac06b 100644
--- a/src/kernels/EnergyExcitation.C
+++ b/src/kernels/EnergyExcitation.C
@@ -33,7 +33,11 @@ EnergyExcitation::computeQpResidual()
 Real
 EnergyExcitation::computeQpJacobian()
 {
-  return 0;
+  return _test[_i][_qp] * _energy_exchange * _N_gas * _electron_density[_qp] *
+         (5e-15 * 0.74 * std::pow((2 / 3), 0.74) * std::pow(_u[_qp], -0.26) *
+              std::exp(-11.56 * 3 / (2 * _u[_qp])) +
+          5e-15 * std::pow((2 * _u[_qp] / 3), 0.74) * (11.56 * 3 / (2 * _u[_qp] * _u[_qp])) *
+              std::exp(-11.56 * 3 / (2 * _u[_qp])));
 }
 
 Real
diff --git a/src/kernels/EnergyIonization.C b/src/kernels/EnergyIonization.C
index 9992457..7f39572 100644
--- a/src/kernels/EnergyIonization.C
+++ b/src/kernels/EnergyIonization.C
@@ -33,7 +33,11 @@ EnergyIonization::computeQpResidual()
 Real
 EnergyIonization::computeQpJacobian()
 {
-  return 0;
+  return _test[_i][_qp] * _energy_exchange * _N_gas * _electron_density[_qp] *
+         (2.34e-14 * 0.59 * std::pow((2 / 3), 0.59) * std::pow(_u[_qp], -0.41) *
+              std::exp(-17.44 * 3 / (2 * _u[_qp])) +
+          2.34e-14 * std::pow((2 * _u[_qp] / 3), 0.59) * (17.44 * 3 / (2 * _u[_qp] * _u[_qp])) *
+              std::exp(-17.44 * 3 / (2 * _u[_qp])));
 }
 
 Real
diff --git a/src/kernels/JouleHeating.C b/src/kernels/JouleHeating.C
index 78e2890..8b18f5b 100644
--- a/src/kernels/JouleHeating.C
+++ b/src/kernels/JouleHeating.C
@@ -53,9 +53,8 @@ JouleHeating::computeQpOffDiagJacobian(unsigned int jvar)
 
   else if (jvar == _potential_id)
     return -_test[_i][_qp] * _e *
-           (-_diffusivity * _grad_electron_density[_qp] +
-            _mobility * _electron_density[_qp] * _grad_phi[_j][_qp]) *
-           _grad_phi[_j][_qp];
+           (-_diffusivity * _grad_electron_density[_qp] * _grad_phi[_j][_qp] +
+            2 * _mobility * _electron_density[_qp] * _grad_potential[_qp] * _grad_phi[_j][_qp]);
 
   else
     return 0;