Update BOUT++ to latest version and fix sheath_ion_heat_flux

src/evolve_density.cxx

@@ -1,16 +1,17 @@
 
 #include <bout/constants.hxx>
-#include <bout/fv_ops.hxx>
-#include <bout/field_factory.hxx>
-#include <bout/output_bout_types.hxx>
 #include <bout/derivs.hxx>
 #include <bout/difops.hxx>
+#include <bout/field_factory.hxx>
+#include <bout/fv_ops.hxx>
 #include <bout/initialprofiles.hxx>
+#include <bout/output_bout_types.hxx>
+#include <bout/solver.hxx>
 
 #include "../include/div_ops.hxx"
 #include "../include/evolve_density.hxx"
-#include "../include/hermes_utils.hxx"
 #include "../include/hermes_build_config.hxx"
+#include "../include/hermes_utils.hxx"
 
 using bout::globals::mesh;
 
@@ -28,9 +29,11 @@ EvolveDensity::EvolveDensity(std::string name, Options& alloptions, Solver* solv
 
   density_floor = options["density_floor"].doc("Minimum density floor").withDefault(1e-7);
 
-  BoutReal temperature_floor = options["temperature_floor"].doc("Low temperature scale for low_T_diffuse_perp")
-    .withDefault<BoutReal>(0.1) / get<BoutReal>(alloptions["units"]["eV"]);
-
+  BoutReal temperature_floor = options["temperature_floor"]
+                                   .doc("Low temperature scale for low_T_diffuse_perp")
+                                   .withDefault<BoutReal>(0.1)
+                               / get<BoutReal>(alloptions["units"]["eV"]);
+
   low_n_diffuse = options["low_n_diffuse"]
                       .doc("Parallel diffusion at low density")
                       .withDefault<bool>(false);
@@ -84,32 +87,33 @@ EvolveDensity::EvolveDensity(std::string name, Options& alloptions, Solver* solv
 
   auto& n_options = alloptions[std::string("N") + name];
   source_time_dependent = n_options["source_time_dependent"]
-    .doc("Use a time-dependent source?")
-    .withDefault<bool>(false);
+                              .doc("Use a time-dependent source?")
+                              .withDefault<bool>(false);
 
   source_only_in_core = n_options["source_only_in_core"]
-    .doc("Zero the source outside the closed field-line region?")
-    .withDefault<bool>(false);
+                            .doc("Zero the source outside the closed field-line region?")
+                            .withDefault<bool>(false);
 
   source_normalisation = Nnorm * Omega_ci;
-  time_normalisation = 1./Omega_ci;
+  time_normalisation = 1. / Omega_ci;
 
   // Try to read the density source from the mesh
   // Units of particles per cubic meter per second
   source = 0.0;
   mesh->get(source, std::string("N") + name + "_src");
   // Allow the user to override the source from input file
   source = n_options["source"]
-    .doc("Source term in ddt(N" + name + std::string("). Units [m^-3/s]"))
-    .withDefault(source)
-    / source_normalisation;
+               .doc("Source term in ddt(N" + name + std::string("). Units [m^-3/s]"))
+               .withDefault(source)
+           / source_normalisation;
 
   // If time dependent, parse the function with respect to time from the input file
   if (source_time_dependent) {
     auto str = n_options["source_prefactor"]
-      .doc("Time-dependent function of multiplier on ddt(N" + name + std::string(") source."))
-      .as<std::string>();
-      source_prefactor_function = FieldFactory::get()->parse(str, &n_options);
+                   .doc("Time-dependent function of multiplier on ddt(N" + name
+                        + std::string(") source."))
+                   .as<std::string>();
+    source_prefactor_function = FieldFactory::get()->parse(str, &n_options);
   }
 
   // Putting source at first X index would put it in both PFR in core, this ensures only core
@@ -126,9 +130,10 @@ EvolveDensity::EvolveDensity(std::string name, Options& alloptions, Solver* solv
     }
   }
 
-  neumann_boundary_average_z = alloptions[std::string("N") + name]["neumann_boundary_average_z"]
-    .doc("Apply neumann boundary with Z average?")
-    .withDefault<bool>(false);
+  neumann_boundary_average_z =
+      alloptions[std::string("N") + name]["neumann_boundary_average_z"]
+          .doc("Apply neumann boundary with Z average?")
+          .withDefault<bool>(false);
 
   std::vector<std::string> outputs = {"AA", "density", "density_source"};
   if (charge != 0.) {
@@ -213,7 +218,9 @@ void EvolveDensity::transform_impl(GuardedOptions& state) {
   if (source_time_dependent) {
     // Evaluate the source_prefactor function at the current time in seconds and scale source with it
     BoutReal time = get<BoutReal>(state["time"]);
-    BoutReal source_prefactor = source_prefactor_function ->generate(bout::generator::Context().set("x",0,"y",0,"z",0,"t",time*time_normalisation));
+    BoutReal source_prefactor =
+        source_prefactor_function->generate(bout::generator::Context().set(
+            "x", 0, "y", 0, "z", 0, "t", time * time_normalisation));
     final_source = source * source_prefactor;
   } else {
     final_source = source;
@@ -230,7 +237,8 @@ void EvolveDensity::finally(const Options& state) {
   // but retain densities which fall below zero
   N.setBoundaryTo(get<Field3D>(species["density"]));
 
-  if ((fabs(charge) > 1e-5) and state.isSection("fields") and state["fields"].isSet("phi")) {
+  if ((fabs(charge) > 1e-5) and state.isSection("fields")
+      and state["fields"].isSet("phi")) {
     // Electrostatic potential set and species is charged -> include ExB flow
 
     Field3D phi = get<Field3D>(state["fields"]["phi"]);
@@ -373,25 +381,26 @@ void EvolveDensity::outputVars(Options& state) {
     auto rho_s0 = get<BoutReal>(state["rho_s0"]);
 
     if (flow_xlow.isAllocated()) {
-      set_with_attrs(state[fmt::format("pf{}_tot_xlow", name)], flow_xlow,
-                   {{"time_dimension", "t"},
-                    {"units", "s^-1"},
-                    {"conversion", rho_s0 * SQ(rho_s0) * Nnorm * Omega_ci},
-                    {"standard_name", "particle flow"},
-                    {"long_name", name + " particle flow in X. Note: May be incomplete."},
-                    {"species", name},
-                    {"source", "evolve_density"}});
+      set_with_attrs(
+          state[fmt::format("pf{}_tot_xlow", name)], flow_xlow,
+          {{"time_dimension", "t"},
+           {"units", "s^-1"},
+           {"conversion", rho_s0 * SQ(rho_s0) * Nnorm * Omega_ci},
+           {"standard_name", "particle flow"},
+           {"long_name", name + " particle flow in X. Note: May be incomplete."},
+           {"species", name},
+           {"source", "evolve_density"}});
     }
     if (flow_ylow.isAllocated()) {
-      set_with_attrs(state[fmt::format("pf{}_tot_ylow", name)], flow_ylow,
-                   {{"time_dimension", "t"},
-                    {"units", "s^-1"},
-                    {"conversion", rho_s0 * SQ(rho_s0) * Nnorm * Omega_ci},
-                    {"standard_name", "particle flow"},
-                    {"long_name", name + " particle flow in Y. Note: May be incomplete."},
-                    {"species", name},
-                    {"source", "evolve_density"}});
+      set_with_attrs(
+          state[fmt::format("pf{}_tot_ylow", name)], flow_ylow,
+          {{"time_dimension", "t"},
+           {"units", "s^-1"},
+           {"conversion", rho_s0 * SQ(rho_s0) * Nnorm * Omega_ci},
+           {"standard_name", "particle flow"},
+           {"long_name", name + " particle flow in Y. Note: May be incomplete."},
+           {"species", name},
+           {"source", "evolve_density"}});
     }
   }
 }
-

src/evolve_energy.cxx

@@ -6,11 +6,12 @@
 #include <bout/initialprofiles.hxx>
 #include <bout/invert_pardiv.hxx>
 #include <bout/output_bout_types.hxx>
+#include <bout/solver.hxx>
 
 #include "../include/div_ops.hxx"
 #include "../include/evolve_energy.hxx"
-#include "../include/hermes_utils.hxx"
 #include "../include/hermes_build_config.hxx"
+#include "../include/hermes_utils.hxx"
 
 #include <algorithm>
 
@@ -226,8 +227,8 @@ void EvolveEnergy::finally(const Options& state) {
 
   Field3D Pfloor = P;
 
-  if (species.isSet("charge") and (fabs(get<BoutReal>(species["charge"])) > 1e-5) and
-      state.isSection("fields") and state["fields"].isSet("phi")) {
+  if (species.isSet("charge") and (fabs(get<BoutReal>(species["charge"])) > 1e-5)
+      and state.isSection("fields") and state["fields"].isSet("phi")) {
     // Electrostatic potential set -> include ExB flow
 
     Field3D phi = get<Field3D>(state["fields"]["phi"]);
@@ -279,7 +280,8 @@ void EvolveEnergy::finally(const Options& state) {
   }
 #if CHECKLEVEL >= 1
   if (species.isSet("pressure_source")) {
-    throw BoutException("Components must evolve `energy_source` rather then `pressure_source`");
+    throw BoutException(
+        "Components must evolve `energy_source` rather then `pressure_source`");
   }
 #endif
   if (species.isSet("momentum_source")) {
@@ -384,24 +386,26 @@ void EvolveEnergy::outputVars(Options& state) {
                     {"source", "evolve_energy"}});
 
     if (flow_xlow.isAllocated()) {
-      set_with_attrs(state[fmt::format("ef{}_tot_xlow", name)], flow_xlow,
-                   {{"time_dimension", "t"},
-                    {"units", "W"},
-                    {"conversion", rho_s0 * SQ(rho_s0) * Pnorm * Omega_ci},
-                    {"standard_name", "power"},
-                    {"long_name", name + " power through X cell face. Note: May be incomplete."},
-                    {"species", name},
-                    {"source", "evolve_energy"}});
+      set_with_attrs(
+          state[fmt::format("ef{}_tot_xlow", name)], flow_xlow,
+          {{"time_dimension", "t"},
+           {"units", "W"},
+           {"conversion", rho_s0 * SQ(rho_s0) * Pnorm * Omega_ci},
+           {"standard_name", "power"},
+           {"long_name", name + " power through X cell face. Note: May be incomplete."},
+           {"species", name},
+           {"source", "evolve_energy"}});
     }
     if (flow_ylow.isAllocated()) {
-      set_with_attrs(state[fmt::format("ef{}_tot_ylow", name)], flow_ylow,
-                   {{"time_dimension", "t"},
-                    {"units", "W"},
-                    {"conversion", rho_s0 * SQ(rho_s0) * Pnorm * Omega_ci},
-                    {"standard_name", "power"},
-                    {"long_name", name + " power through Y cell face. Note: May be incomplete."},
-                    {"species", name},
-                    {"source", "evolve_energy"}});
+      set_with_attrs(
+          state[fmt::format("ef{}_tot_ylow", name)], flow_ylow,
+          {{"time_dimension", "t"},
+           {"units", "W"},
+           {"conversion", rho_s0 * SQ(rho_s0) * Pnorm * Omega_ci},
+           {"standard_name", "power"},
+           {"long_name", name + " power through Y cell face. Note: May be incomplete."},
+           {"species", name},
+           {"source", "evolve_energy"}});
     }
   }
 }