Merge branch 'neutral-advection' into fix-2d-neutral-eqns

bendudson · May 28, 2024 · 6ede2cb · 6ede2cb
2 parents ee46734 + c095e43
commit 6ede2cb
Show file tree

Hide file tree

Showing 2 changed files with 31 additions and 4 deletions.
diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx
@@ -75,6 +75,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
@@ -391,17 +391,23 @@ void NeutralMixed::finally(const Options& state) {
     + Div_a_Grad_perp_upwind_flows((5. / 3) * DnnPn, logPnlim,          // Perpendicular advection
                                    energy_flow_xlow, energy_flow_ylow)  
      ;
-  energy_flow_xlow *= 3/2; // Note: Should this be 5/2?
-  energy_flow_ylow *= 3/2;
+  // The factor here is likely 5/2 as we're advecting internal energy and pressure.
+  // Doing this still leaves a heat imbalance factor of 0.11 in the cells, but better than 0.33 with 3/2.
+  energy_flow_xlow *= 5/2; 
+  energy_flow_ylow *= 5/2;
 
 
 
   if (neutral_conduction) {
-    ddt(Pn) += 
-      (2. / 3) * Div_a_Grad_perp_upwind(kappa_n, Tn)                      // Perpendicular conduction
+    ddt(Pn) += (2. / 3) Div_a_Grad_perp_upwind_flows(kappa_n, Tn,
+                        conduction_flow_xlow, conduction_flow_ylow)    // Perpendicular conduction
       + FV::Div_par_K_Grad_par(kappa_n, Tn)                               // Parallel conduction
       ;
   }
+
+    // The factor here is likely 3/2 as this is pure energy flow, but needs checking.
+  conduction_flow_xlow *= 3/2;
+  conduction_flow_ylow *= 3/2;
 
   Sp = pressure_source;
   if (localstate.isSet("energy_source")) {
@@ -661,6 +667,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"}});
+    }
   }
 }