Merge pull request #385 from boutproject/more-diagnostics

bendudson · web-flow · commit a35e2f5e3ec7 · 2025-08-27T21:06:15.000-07:00
More diagnostics
diff --git a/include/amjuel_reaction.hxx b/include/amjuel_reaction.hxx
@@ -146,6 +146,13 @@ protected:
       add(electron["density_source"], (to_charge - from_charge) * reaction_rate);
       if (electron.isSet("velocity")) {
         // Transfer of electron kinetic to thermal energy due to density source
+        // For ionisation:
+        // Electrons with zero average velocity are created, diluting the kinetic energy.
+        // Total energy conservation requires a corresponding internal energy source.
+        //
+        // For recombination:
+        // Electrons with some velocity are incorporated into a neutral species with their kinetic
+        // energy converted to an internal energy source of that species.
         auto Ve = get<Field3D>(electron["velocity"]);
         auto Ae = get<BoutReal>(electron["AA"]);
         add(electron["energy_source"], 0.5 * Ae * (to_charge - from_charge) * reaction_rate * SQ(Ve));
@@ -182,6 +189,11 @@ protected:
     // d/dt(3/2 p_2) = - m R v_1 v_2 + (1/2) m R v_1^2 + (1/2) m R v_2^2
     //               = (1/2) m R (v_1 - v_2)^2
     //
+    // This term accounts for the broadening of the species velocity
+    // distribution when combining the two underlying distributions. 
+    // The greater the difference in velocities of the underlying species,
+    // the wider the resultant distribution, which corresponds to an
+    // increase in temperature and therefore internal energy.
 
     add(to_ion["energy_source"], 0.5 * AA * reaction_rate * SQ(V1 - V2));
 
diff --git a/include/collisions.hxx b/include/collisions.hxx
@@ -61,7 +61,7 @@ private:
 
   /// Calculated collision rates saved for post-processing and use by other components
   /// Saved in options, the BOUT++ dictionary-like object
-  Options collision_rates;
+  Options collision_rates, energy_channels, friction_energy_channels, momentum_channels;
 
   /// Save more diagnostics?
   bool diagnose;
diff --git a/include/evolve_pressure.hxx b/include/evolve_pressure.hxx
@@ -116,11 +116,13 @@ private:
   bool source_time_dependent; ///< Is the input source time dependent?
   Field3D flow_xlow, flow_ylow; ///< Energy flow diagnostics
   Field3D flow_ylow_conduction; ///< Conduction energy flow diagnostics
-  Field3D flow_ylow_kinetic;    ///< Parallel flow of kinetic energy
+  Field3D flow_ylow_advection;    ///< Advection energy flow diagnostics
+  Field3D flow_ylow_viscous_heating; ///< Flow of kinetic energy due to numerical viscosity
 
   bool numerical_viscous_heating; ///< Include heating due to numerical viscosity?
   bool fix_momentum_boundary_flux; ///< Fix momentum flux to boundary condition?
   Field3D Sp_nvh; ///< Pressure source due to artificial viscosity
+  Field3D E_PdivV, E_VgradP; ///< Diagnostic energy source terms for p*Div(V) and V*Grad(P)
 };
 
 namespace {
diff --git a/src/collisions.cxx b/src/collisions.cxx
@@ -104,6 +104,10 @@ void Collisions::collide(Options& species1, Options& species2, const Field3D& nu
 
       add(species1["momentum_source"], F12);
       subtract(species2["momentum_source"], F12);
+      
+      // Diagnostics
+      set(momentum_channels[species1.name()][species2.name()], F12); 
+      set(momentum_channels[species2.name()][species1.name()], -F12);
 
       if (frictional_heating) {
         // Heating due to friction and energy transfer
@@ -127,8 +131,15 @@ void Collisions::collide(Options& species1, Options& species2, const Field3D& nu
         //  1) This term is always positive: Collisions don't lead to cooling
         //  2) In the limit that m_2 << m_1 (e.g. electron-ion collisions),
         //     the lighter species is heated more than the heavy species.
-        add(species1["energy_source"], (A2 / (A1 + A2)) * (velocity2 - velocity1) * F12);
-        add(species2["energy_source"], (A1 / (A1 + A2)) * (velocity2 - velocity1) * F12);
+        Field3D species1_source = (A2 / (A1 + A2)) * (velocity2 - velocity1) * F12;
+        Field3D species2_source = (A1 / (A1 + A2)) * (velocity2 - velocity1) * F12;
+
+        add(species1["energy_source"], species1_source);
+        add(species2["energy_source"], species2_source);
+
+        // Diagnostics
+        set(friction_energy_channels[species1.name()][species2.name()], species1_source);
+        set(friction_energy_channels[species2.name()][species1.name()], species2_source);
       }
     }
 
@@ -144,6 +155,10 @@ void Collisions::collide(Options& species1, Options& species2, const Field3D& nu
 
       add(species1["energy_source"], Q12);
       subtract(species2["energy_source"], Q12);
+
+      // Diagnostics
+      set(energy_channels[species1.name()][species2.name()], Q12);
+      set(energy_channels[species2.name()][species1.name()], -Q12);
     }
   }
 }
@@ -488,6 +503,10 @@ void Collisions::outputVars(Options& state) {
 
   // Normalisations
   auto Omega_ci = get<BoutReal>(state["Omega_ci"]);
+  auto Nnorm = get<BoutReal>(state["Nnorm"]);
+  auto Tnorm = get<BoutReal>(state["Tnorm"]);
+  BoutReal Pnorm = SI::qe * Tnorm * Nnorm; // Pressure normalisation
+  auto Cs0 = get<BoutReal>(state["Cs0"]);
 
   /// Iterate through the first species in each collision pair
   const std::map<std::string, Options>& level1 = collision_rates.getChildren();
@@ -498,17 +517,65 @@ void Collisions::outputVars(Options& state) {
     const std::map<std::string, Options>& level2 = section.getChildren();
     for (auto s2 = std::begin(level2); s2 != std::end(level2); ++s2) {
 
-      std::string name = s1->first + s2->first;
+      std::string A = s1->first;
+      std::string B = s2->first;
+      std::string AB = A + B;
 
-      set_with_attrs(state[std::string("K") + name + std::string("_coll")],
+      // Collision frequencies
+      set_with_attrs(state[std::string("K") + AB + std::string("_coll")],
                      getNonFinal<Field3D>(section[s2->first]),
                      {{"time_dimension", "t"},
                       {"units", "s-1"},
                       {"conversion", Omega_ci},
-                      {"standard_name", "collision frequency"},
-                      {"long_name", name + " collision frequency"},
-                      {"species", name},
+                      {"standard_name", AB + "collision frequency"},
+                      {"long_name", AB + " collision frequency"},
+                      {"species", A},
                       {"source", "collisions"}});
+
+      // Collisional energy transfer channels (i.e. thermal equilibration)
+      if ((energy_channels.isSection(A)) and (energy_channels[A].isSet(B))) {
+
+        set_with_attrs(state[std::string("E") + AB + std::string("_coll")],
+                     getNonFinal<Field3D>(energy_channels[A][B]),
+                     {{"time_dimension", "t"},
+                      {"units", "W / m^3"},
+                      {"conversion", Pnorm * Omega_ci},
+                      {"standard_name", AB + "collisional energy transfer source"},
+                      {"long_name", AB + "collisional energy transfer source"},
+                      {"species", A},
+                      {"source", "collisions"}});
+      }
+
+      // Frictional energy sources (both species heat through friction)
+      if ((friction_energy_channels.isSection(A)) and (friction_energy_channels[A].isSet(B))) {
+
+        set_with_attrs(state[std::string("E") + AB + std::string("_coll_friction")],
+                     getNonFinal<Field3D>(friction_energy_channels[A][B]),
+                     {{"time_dimension", "t"},
+                      {"units", "W / m^3"},
+                      {"conversion", Pnorm * Omega_ci},
+                      {"standard_name", AB + "frictional energy source"},
+                      {"long_name", AB + "frictional energy source"},
+                      {"species", A},
+                      {"source", "collisions"}});
+      }
+
+      // Momentum exchange channels
+      if ((momentum_channels.isSection(A)) and (momentum_channels[A].isSet(B))) {
+
+        set_with_attrs(state[std::string("F") + AB + std::string("_coll")],
+                     getNonFinal<Field3D>(friction_energy_channels[A][B]),
+                     {{"time_dimension", "t"},
+                      {"units", "kg m^-2 s^-2"},
+                      {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci},
+                      {"standard_name", AB + "collisional momentum transfer"},
+                      {"long_name", AB + "collisional momentum transfer"},
+                      {"species", A},
+                      {"source", "collisions"}});
+      }
+
     }
   }
+
+
 }
diff --git a/src/evolve_pressure.cxx b/src/evolve_pressure.cxx
@@ -289,21 +289,24 @@ void EvolvePressure::finally(const Options& state) {
 
     if (p_div_v) {
       // Use the P * Div(V) form
-      ddt(P) -= FV::Div_par_mod<hermes::Limiter>(P, V, fastest_wave, flow_ylow);
+      ddt(P) -= FV::Div_par_mod<hermes::Limiter>(P, V, fastest_wave, flow_ylow_advection);
 
       // Work done. This balances energetically a term in the momentum equation
-      ddt(P) -= (2. / 3) * Pfloor * Div_par(V);
+      E_PdivV = -Pfloor * Div_par(V);
+      ddt(P) += (2. / 3) * E_PdivV;
 
     } else {
       // Use V * Grad(P) form
       // Note: A mixed form has been tried (on 1D neon example)
       //       -(4/3)*FV::Div_par(P,V) + (1/3)*(V * Grad_par(P) - P * Div_par(V))
       //       Caused heating of charged species near sheath like p_div_v
-      ddt(P) -= (5. / 3) * FV::Div_par_mod<hermes::Limiter>(P, V, fastest_wave, flow_ylow);
+      ddt(P) -= (5. / 3) * FV::Div_par_mod<hermes::Limiter>(P, V, fastest_wave, flow_ylow_advection);
 
-      ddt(P) += (2. / 3) * V * Grad_par(P);
+      E_VgradP =  V * Grad_par(P);
+      ddt(P) += (2. / 3) * E_VgradP;
     }
-    flow_ylow *= 5. / 2; // Energy flow
+    flow_ylow_advection *= 5. / 2; // Energy flow
+    flow_ylow = flow_ylow_advection;
 
     if (state.isSection("fields") and state["fields"].isSet("Apar_flutter")) {
       // Magnetic flutter term
@@ -313,12 +316,12 @@ void EvolvePressure::finally(const Options& state) {
     }
 
     if (numerical_viscous_heating || diagnose) {
-      // Viscous heating coming from numerical viscosity
+      // Viscous heating coming from numerical viscosity from the Lax flux
       Field3D Nlim = softFloor(N, density_floor);
       const BoutReal AA = get<BoutReal>(species["AA"]); // Atomic mass
-      Sp_nvh = (2. / 3) * AA * FV::Div_par_fvv_heating(Nlim, V, fastest_wave, flow_ylow_kinetic, fix_momentum_boundary_flux);
-      flow_ylow_kinetic *= AA;
-      flow_ylow += flow_ylow_kinetic;
+      Sp_nvh = (2. / 3) * AA * FV::Div_par_fvv_heating(Nlim, V, fastest_wave, flow_ylow_viscous_heating, fix_momentum_boundary_flux);
+      flow_ylow_viscous_heating *= AA;
+      flow_ylow += flow_ylow_viscous_heating;
       if (numerical_viscous_heating) {
         ddt(P) += Sp_nvh;
       }
@@ -652,6 +655,30 @@ void EvolvePressure::outputVars(Options& state) {
                     {"species", name},
                     {"source", "evolve_pressure"}});
 
+    if (p_div_v) {
+
+      set_with_attrs(state["E" + name + "_PdivV"], E_PdivV,
+                   {{"time_dimension", "t"},
+                    {"units", "W / m^-3"},
+                    {"conversion", Pnorm * Omega_ci},
+                    {"standard_name", "energy source"},
+                    {"long_name", name + " energy source due to pressure gradient"},
+                    {"species", name},
+                    {"source", "evolve_pressure"}});
+    } else {
+
+      set_with_attrs(state["E" + name + "_VgradP"], E_VgradP,
+                   {{"time_dimension", "t"},
+                    {"units", "W / m^-3"},
+                    {"conversion", Pnorm * Omega_ci},
+                    {"standard_name", "energy source"},
+                    {"long_name", name + " energy source due to pressure gradient"},
+                    {"species", name},
+                    {"source", "evolve_pressure"}});
+
+    }
+                  
+
     if (flow_xlow.isAllocated()) {
       set_with_attrs(state[fmt::format("ef{}_tot_xlow", name)], flow_xlow,
                    {{"time_dimension", "t"},
@@ -679,22 +706,34 @@ void EvolvePressure::outputVars(Options& state) {
                     {"units", "W"},
                     {"conversion", rho_s0 * SQ(rho_s0) * Pnorm * Omega_ci},
                     {"standard_name", "power"},
-                    {"long_name", name + " conduction through Y cell face. Note: May be incomplete."},
+                    {"long_name", name + " conduction through Y cell face."},
                     {"species", name},
                     {"source", "evolve_pressure"}});
     }
 
-    if (flow_ylow_kinetic.isAllocated()) {
-      set_with_attrs(state[fmt::format("ef{}_kin_ylow", name)], flow_ylow_kinetic,
+    if (flow_ylow_advection.isAllocated()) {
+      set_with_attrs(state[fmt::format("ef{}_adv_ylow", name)], flow_ylow_advection,
                    {{"time_dimension", "t"},
                     {"units", "W"},
                     {"conversion", rho_s0 * SQ(rho_s0) * Pnorm * Omega_ci},
                     {"standard_name", "power"},
-                    {"long_name", name + " kinetic energy flow through Y cell face. Note: May be incomplete."},
+                    {"long_name", name + " advected energy flow through Y cell face."},
                     {"species", name},
                     {"source", "evolve_pressure"}});
     }
 
+    if (flow_ylow_viscous_heating.isAllocated()) {
+      set_with_attrs(state[fmt::format("ef{}_visc_heat_ylow", name)], flow_ylow_viscous_heating,
+                   {{"time_dimension", "t"},
+                    {"units", "W"},
+                    {"conversion", rho_s0 * SQ(rho_s0) * Pnorm * Omega_ci},
+                    {"standard_name", "power"},
+                    {"long_name", name + " energy flow due to Lax flux numerical viscosity through Y cell face"},
+                    {"species", name},
+                    {"source", "evolve_pressure"}});
+    }
+
+
     if (numerical_viscous_heating) {
       set_with_attrs(state[std::string("E") + name + std::string("_nvh")], Sp_nvh * 3/.2,
                    {{"time_dimension", "t"},
