add excitation to the radiation sink term

Author: hahahasan

src/full-velocity.cxx

@@ -324,10 +324,10 @@ void FullVelocity::update(const Field3D &Ne, const Field3D &Te,
   /////////////////////////////////////////////////////
   // Atomic processes
 
-  Field3D Riz, Rrc, Rcx;
-  neutral_rates(Ne, Te, Ti, Vi, Nn, Tn, Vnpar, S, F, Qi, Rp, Riz, Rrc, Rcx);
+  Field3D Riz, Rrc, Rcx, Rex;
+  neutral_rates(Ne, Te, Ti, Vi, Nn, Tn, Vnpar, S, F, Qi, Rp, Riz, Rrc, Rcx, Rex);
 
-  Fperp = Rrc + Rcx; // Friction for vorticity
+  Fperp = Rrc + Rcx + Rex; // Friction for vorticity
 
   // Loss of momentum in the X and Z directions
   ddt(Vn2D).x -= (DC(Rcx) + DC(Riz)) * Vn2D.x / Nn2D_floor;

src/mixed.cxx

@@ -168,13 +168,13 @@ void NeutralMixed::update(const Field3D &Ne, const Field3D &Te,
   // Atomic processes
   TRACE("Atomic processes");
 
-  Field3D Riz, Rrc, Rcx;
-  neutral_rates(Ne, Te, Ti, Vi, Nn, Tn, Vn, S, F, Qi, Rp, Riz, Rrc, Rcx);
+  Field3D Riz, Rrc, Rcx, Rex;
+  neutral_rates(Ne, Te, Ti, Vi, Nn, Tn, Vn, S, F, Qi, Rp, Riz, Rrc, Rcx, Rex);
 
   // Neutral cross-field diffusion coefficient
   BoutReal neutral_lmax = 0.1 / Lnorm;
   Field3D Rnn = Nn * sqrt(Tn) / neutral_lmax; // Neutral-neutral collisions
-  Dnn = Pnlim / (Riz + Rcx + Rnn);
+  Dnn = Pnlim / (Riz + Rcx + Rnn + Rex);
 
   mesh->communicate(Dnn);
   Dnn.clearParallelSlices();
@@ -254,7 +254,7 @@ void NeutralMixed::update(const Field3D &Ne, const Field3D &Te,
       + FV::Div_par_K_Grad_par((2. / 5) * DnnNn, Vn) // Parallel viscosity
       ;
 
-  Fperp = Rcx + Rrc;
+  Fperp = Rcx + Rrc + Rex;
 
   /////////////////////////////////////////////////////
   // Neutral pressure

src/neutral-model.cxx

@@ -48,7 +48,7 @@ void NeutralModel::neutral_rates(
     const Field3D &Vi, // Plasma quantities
     const Field3D &Nn, const Field3D &Tn, const Field3D &Vnpar, // Neutral gas
     Field3D &S, Field3D &F, Field3D &Qi, Field3D &R, // Transfer rates
-    Field3D &Riz, Field3D &Rrc, Field3D &Rcx) {      // Rates
+    Field3D &Riz, Field3D &Rrc, Field3D &Rcx, Field3D &Rex) {      // Rates
 
   // Allocate output fields
   S = 0.0;
@@ -184,6 +184,27 @@ void NeutralModel::neutral_rates(
         BoutReal R_iz_R =
             Ne_R * Nn_R * hydrogen.ionisation(Te_R * Tnorm) * Nnorm / Fnorm;
 
+        // Excitation
+        if (excitation) {
+          /////////////////////////////////////////////////////////
+          // Electron-neutral excitation
+          // Note: Rates need checking
+          // Currently assuming that quantity calculated is in [eV m^3/s]
+
+          BoutReal R_ex_L = Ne_L * Nn_L *
+                            hydrogen.excitation(Te_L * Tnorm) * Nnorm /
+                            Fnorm / Tnorm;
+          BoutReal R_ex_C = Ne_C * Nn_C *
+                            hydrogen.excitation(Te_C * Tnorm) * Nnorm /
+                            Fnorm / Tnorm;
+          BoutReal R_ex_R = Ne_R * Nn_R *
+                            hydrogen.excitation(Te_R * Tnorm) * Nnorm /
+                            Fnorm / Tnorm;
+
+          Rex(i, j, k) = (J_L * R_ex_L + 4. * J_C * R_ex_C + J_R * R_ex_R) /
+                          (6. * J_C);
+        }
+
         // Neutral sink, plasma source
         S(i, j, k) -=
             (J_L * R_iz_L + 4. * J_C * R_iz_C + J_R * R_iz_R) / (6. * J_C);
@@ -200,6 +221,8 @@ void NeutralModel::neutral_rates(
                        (6. * J_C);
 
         // Ionisation and electron excitation energy
+
+
         R(i, j, k) += (Eionize / Tnorm) *
                       (J_L * R_iz_L + 4. * J_C * R_iz_C + J_R * R_iz_R) /
                       (6. * J_C);