Merge pull request #167 from AetherModel/cet

FEAT: Electron Temperature

Author: aaronjridley

include/inputs.h

@@ -436,7 +436,50 @@ class Inputs {
   bool get_O_cooling();
 
   /**********************************************************************
-      \brief returns settings["
+     \brief returns settings["Sources"]["Ions"]["IncludePhotoElectronHeating"]
+     \param none
+   **/
+  bool get_do_photoelectron_heating();
+
+  /**********************************************************************
+     \brief returns settings["Sources"]["Ions"]["IncludeIonizationHeating"]
+     \param none
+   **/
+  bool get_do_ionization_heating();
+
+  /**********************************************************************
+     \brief returns settings["Sources"]["Ions"]["IncludeElectronIonCollisionalHeating"]
+     \param none
+   **/
+  bool get_do_electron_ion_collisional_heating();
+
+  /**********************************************************************
+     \brief returns settings["Sources"]["Ions"]["IncludeElectronNeutralElasticCollisionalHeating"]
+     \param none
+   **/
+  bool get_do_electron_neutral_elastic_collisional_heating();
+
+  /**********************************************************************
+     \brief returns settings["Sources"]["Ions"]["IncludeElectronNeutralInelasticCollisionalHeating"]
+     \param none
+   **/
+  bool get_do_electron_neutral_inelastic_collisional_heating();
+
+
+  /**********************************************************************
+     \brief returns settings["Sources"]["Ions"]["IncludeThermoelectricHeating"]
+     \param none
+   **/
+  bool get_do_thermoelectric_heating();
+
+  /**********************************************************************
+     \brief returns settings["
+     \param 
+   **/
+  bool get_use_centripetal();
+  
+  /**********************************************************************
+     \brief returns settings["
      \param 
    **/
   bool get_use_coriolis();

include/ions.h

@@ -172,7 +172,7 @@ class Ions {
   /// Specific heat (constant volume):
   arma_cube Cv_scgc;
 
-  /// Head Conduction (bulk):
+  /// Heat Conduction (bulk):
   arma_cube lambda;
 
   // Electrodynamics:
@@ -348,8 +348,75 @@ class Ions {
      \brief Calculates the electron temperature on the given grid
      \param neutrals these are needed for the collision terms
      \param grid this is the grid to solve the equation on
+     \param time the time class to know dt
    **/
-  void calc_electron_temperature(Neutrals neutrals, Grid grid);
+  void calc_electron_temperature(Neutrals neutrals, Grid grid, Times time);
+
+  /**********************************************************************
+  /// @brief Calculate epsilon
+  /// @details intermediate variable used in photoelectron & ionization heating
+  /// From (Smithro & Solomon, 2008).
+  /// @param neutrals 
+  /// @return epsilon
+  **/
+  arma_cube calc_epsilon(Neutrals &neutrals);
+
+  /**********************************************************************
+    \brief Calculates photoelectron heating
+    \details Based on (Swartz & Nisbet, 1972) & (Smithro & Solomon, 2008)
+            Uses equations 9-12 from (Zhu & Ridley, 2016)
+            https://doi.org/10.1016/j.jastp.2016.01.005
+    \param epsilon 
+    \return Qphe 
+  **/
+  arma_cube calc_photoelectron_heating(arma_cube epsilon);
+
+  /**********************************************************************
+    \brief Calculates auroral heating
+    \details NOTE: in GITM this is solved separately for ion precipitation & auroral 
+        ionization. In Aether these are both in ions.species[iIon].ionization_scgc...
+    \param epsilon 
+    \return Qaurora 
+  **/
+  arma_cube calc_ionization_heating(arma_cube epsilon);
+
+  /**********************************************************************
+    \brief Calculates electron-ion (elastic) collisional heating
+    \details From Schunk and Nagy 2009, and Bei-Chen Zhang and Y. Kamide 2003
+    - This differs slightly from the GITM implementation, which assumes several ion species are present.
+      Instead, here we use each ion species for the sum.
+    - electon-ion collision frequency (from Schunk and Nagy 2009) = 5.45E-5
+    - This is capable of handling BOTH the bulk & individual ion temperatures
+    \return vector<Qeicp, Qeicm, Qeic_v>
+  **/
+  std::vector<arma_cube> calc_electron_ion_collisions();
+
+  /**********************************************************************
+    \brief Calculates electron-neutral elastic collisional heating
+    \details From Schunk and Nagy 2009
+    \param neutrals 
+    \return vector<Qencp, Qencm, Qenc_v>
+  **/
+  std::vector<arma_cube> calc_electron_neutral_elastic_collisions(Neutrals &neutrals);
+
+  /**********************************************************************
+    \brief Calculates the electron-neutral inelastic collisional heating
+    \details From Schunk and Nagy 2009 pages 277, 282.
+    This includes N2, O2 rotation, fine structure, O(1D) exitation & vibration, N2 vibration.
+    See equation 15 from (Zhu, Ridley, Deng, 2016) https://doi.org/10.1016/j.jastp.2016.01.005
+    \param neutrals 
+    \return vector<Qencp, Qencm, Qenc_v>
+  **/
+  std::vector<arma_cube> calc_electron_neutral_inelastic_collisions(Neutrals &neutrals);
+
+  /**********************************************************************
+    \brief Calculate the thermoelectric current (same at all altitudes)
+    \details Use eq. 6 of https://doi.org/10.1016/j.jastp.2016.01.005
+    - Since we do not know e- parallel velocity, the dipole needs to do it this way too.
+    \param grid 
+    \return arma_mat JParaAlt
+  **/
+  arma_mat calc_thermoelectric_current(Grid &grid);
 
   /**********************************************************************
      \brief Check all of the variables for nonfinites, such as nans

include/solvers.h

@@ -40,7 +40,8 @@ arma_vec solver_conduction(
 			   arma_vec dx,
 			   precision_t dt,
 			   int64_t nGCs,
-			   bool return_diff);
+			   bool return_diff = false,
+			   arma_vec source2 = arma_vec());
 
 arma_cube solver_chemistry(arma_cube density,
 			   arma_cube source,

share/run/UA/inputs/defaults.json

@@ -104,7 +104,14 @@
 		"Cent_acc": true },
 	"Neutrals" : {
 	    "NO_cool" : false,
-	    "O_cool": false } },
+	    "O_cool": false },
+	"Ions":{
+		"IncludePhotoElectronHeating": true,
+		"IncludeIonizationHeating": true,
+		"IncludeElectronIonCollisionalHeating":true,
+		"IncludeElectronNeutralElasticCollisionalHeating":true,
+		"IncludeElectronNeutralInelasticCollisionalHeating":true}
+	},
 
     "Seed" : 0,
 

src/advance.cpp

@@ -200,9 +200,9 @@ bool advance(Planets &planet,
     }
 
     ions.calc_ion_temperature(neutrals, gGrid, time);
-    ions.calc_electron_temperature(neutrals, gGrid);
+    // ions.calc_electron_temperature(neutrals, gGrid, time);
     //ionsMag.calc_ion_temperature(neutralsMag, mGrid, time);
-    //ionsMag.calc_electron_temperature(neutralsMag, mGrid);
+    ionsMag.calc_electron_temperature(neutralsMag, mGrid, time);
 
     if (didWork & input.get_check_for_nans())
       didWork = neutrals.check_for_nonfinites("After Vertical Advection");