From 0aa087b0c2b21e6445a0f84400b11c813fff573a Mon Sep 17 00:00:00 2001
From: Mike Kryjak <mikekryjak@gmail.com>
Date: Tue, 28 May 2024 15:21:21 +0100
Subject: [PATCH] Change conduction to upwind and expose fluxes

---
 include/neutral_mixed.hxx |  1 +
 src/neutral_mixed.cxx     | 23 ++++++++++++++++++++++-
 2 files changed, 23 insertions(+), 1 deletion(-)

diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx
index baaff4de..3747878f 100644
--- a/include/neutral_mixed.hxx
+++ b/include/neutral_mixed.hxx
@@ -73,6 +73,7 @@ private:
   Field3D particle_flow_xlow, particle_flow_ylow;
   Field3D momentum_flow_xlow, momentum_flow_ylow;
   Field3D energy_flow_xlow, energy_flow_ylow;
+  Field3D conduction_flow_xlow, conduction_flow_ylow;
 };
 
 namespace {
diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx
index 7d9811eb..3b71815b 100644
--- a/src/neutral_mixed.cxx
+++ b/src/neutral_mixed.cxx
@@ -421,7 +421,8 @@ void NeutralMixed::finally(const Options& state) {
   energy_flow_ylow *= 5/2;
 
   if (neutral_conduction) {
-    ddt(Pn) += FV::Div_a_Grad_perp(DnnNn, Tn)    // Perpendicular conduction
+    ddt(Pn) += Div_a_Grad_perp_upwind_flows(DnnNn, Tn,
+                        conduction_flow_xlow, conduction_flow_ylow)    // Perpendicular conduction
       + FV::Div_par_K_Grad_par(DnnNn, Tn)        // Parallel conduction
       ;
   }
@@ -646,6 +647,26 @@ void NeutralMixed::outputVars(Options& state) {
                     {"species", name},
                     {"source", "evolve_pressure"}});
     }
+    if (conduction_flow_xlow.isAllocated()) {
+      set_with_attrs(state[std::string("ConductionFlow_") + name + std::string("_xlow")],conduction_flow_xlow,
+                   {{"time_dimension", "t"},
+                    {"units", "W"},
+                    {"conversion", rho_s0 * SQ(rho_s0) * Pnorm * Omega_ci},
+                    {"standard_name", "power"},
+                    {"long_name", name + " conducted power through X cell face. Note: May be incomplete."},
+                    {"species", name},
+                    {"source", "evolve_pressure"}});
+    }
+    if (conduction_flow_ylow.isAllocated()) {
+      set_with_attrs(state[std::string("ConductionFlow_") + name + std::string("_ylow")], conduction_flow_ylow,
+                   {{"time_dimension", "t"},
+                    {"units", "W"},
+                    {"conversion", rho_s0 * SQ(rho_s0) * Pnorm * Omega_ci},
+                    {"standard_name", "power"},
+                    {"long_name", name + " conducted power through Y cell face. Note: May be incomplete."},
+                    {"species", name},
+                    {"source", "evolve_pressure"}});
+    }
   }
 }