have outputs for the ionisation, charge-exchange and recombination rates and included the option to reverse the Curl(b/B) vector. Also had some small code tidy ups

Author: hahahasan

include/hermes-2.hxx

@@ -226,6 +226,7 @@ private:
 
   // Curvature, Grad-B drift
   Vector3D Curlb_B; // Curl(b/B)
+  bool revCurlb_B; // Reverse direction of Curl(b/B) vector
 
   // Perturbed parallel gradient operators
   const Field3D Grad_parP(const Field3D &f);

include/mixed.hxx

@@ -33,6 +33,8 @@ private:
   Field3D Pnlim;  // Limited pressure, used to calculate pressure-driven diffusive flows
   Field3D Vn;
   Field3D Dnn;
+
+  Field3D Riz, Rrc, Rcx, Rex;
 
   bool sheath_ydown, sheath_yup;
 

src/hermes-2.cxx

@@ -307,6 +307,10 @@ int Hermes::init(bool restarting) {
   poloidal_flows = optsc["poloidal_flows"]
                        .doc("Include ExB flows in X-Y plane")
                        .withDefault(true);
+
+  revCurlb_B = optsc["revCurlb_B"]
+                    .doc("Reverse the direction of Curl(b/B)")
+                    .withDefault(true);
 
   OPTION(optsc, ion_velocity, true);
 
@@ -522,6 +526,12 @@ int Hermes::init(bool restarting) {
   // Inflowing density carries momentum
   OPTION(optne, density_inflow, false);
 
+  auto& optpe = opt["Pe"];
+  PeSource = optpe["source"].withDefault(Field3D{0.0});
+
+  auto& optpi = opt["Pi"];
+  PiSource = optpi["source"].withDefault(Field3D{0.0});
+
   // radial buffer setup
   if (slab_radial_buffers || radial_buffers) {
     // Need to set the sources in the radial buffer regions to zero
@@ -585,13 +595,13 @@ int Hermes::init(bool restarting) {
   NeSource /= Omega_ci;
   Sn = DC(NeSource);
 
-  auto& optpe = opt["Pe"];
-  PeSource = optpe["source"].withDefault(Field3D{0.0});
+  // auto& optpe = opt["Pe"];
+  // PeSource = optpe["source"].withDefault(Field3D{0.0});
   PeSource /= Omega_ci;
   Spe = DC(PeSource);
 
-  auto& optpi = opt["Pi"];
-  PiSource = optpi["source"].withDefault(Field3D{0.0});
+  // auto& optpi = opt["Pi"];
+  // PiSource = optpi["source"].withDefault(Field3D{0.0});
   PiSource /= Omega_ci;
   Spi = DC(PiSource);
 
@@ -603,7 +613,7 @@ int Hermes::init(bool restarting) {
         for (int y = mesh->ystart; y <= mesh->yend; y++) {
           Sn(x, y) = 0.0;
           Spe(x, y) = 0.0;
-	  Spi(x, y) = 0.0;
+	        Spi(x, y) = 0.0;
         }
       }
     }
@@ -907,6 +917,10 @@ int Hermes::init(bool restarting) {
       }
     }
   }
+
+  if (revCurlb_B) {
+    Curlb_B = -1.0 * Curlb_B
+  }
 
   if (Options::root()["mesh"]["paralleltransform"].as<std::string>() == "shifted") {
     // Check if the gridfile was created for "shiftedmetric" or for "identity" parallel

src/mixed.cxx

@@ -64,6 +64,8 @@ NeutralMixed::NeutralMixed(Solver *solver, Mesh *UNUSED(mesh), Options &options)
   F = 0;
   Qi = 0;
   Rp = 0;
+
+  SAVE_REPEAT4(Riz, Rrc, Rcx, Rex);
 }
 
 void NeutralMixed::update(const Field3D &Ne, const Field3D &Te,
@@ -168,7 +170,6 @@ void NeutralMixed::update(const Field3D &Ne, const Field3D &Te,
   // Atomic processes
   TRACE("Atomic processes");
 
-  Field3D Riz, Rrc, Rcx, Rex;
   neutral_rates(Ne, Te, Ti, Vi, Nn, Tn, Vn, S, F, Qi, Rp, Riz, Rrc, Rcx, Rex);
 
   // Neutral cross-field diffusion coefficient

src/neutral-model.cxx

@@ -61,8 +61,6 @@ void NeutralModel::neutral_rates(
   Rcx = 0.0;
   Rex = 0.0;
 
-  SAVE_REPEAT4(Riz, Rrc, Rcx, Rex);
-
   Coordinates *coord = mesh->getCoordinates();
 
   for (int i = mesh->xstart; i <= mesh->xend; i++)