From 155f10e1d0fce4f8f880c99248302282381806fd Mon Sep 17 00:00:00 2001
From: Mike Kryjak <mikekryjak@gmail.com>
Date: Thu, 1 Aug 2024 09:31:49 +0100
Subject: [PATCH] Add flag to suppress outflow (currently incorrect) - Intent
 is to preserve heat flux but prevent particle flux - Enables test case
 without neutrals

---
 include/sheath_boundary_simple.hxx |  2 ++
 src/sheath_boundary_simple.cxx     | 29 +++++++++++++++++++++++++----
 2 files changed, 27 insertions(+), 4 deletions(-)

diff --git a/include/sheath_boundary_simple.hxx b/include/sheath_boundary_simple.hxx
index e03cb9c2..5f63187c 100644
--- a/include/sheath_boundary_simple.hxx
+++ b/include/sheath_boundary_simple.hxx
@@ -86,6 +86,8 @@ private:
 
   bool always_set_phi; ///< Set phi field?
 
+  bool suppress_outflow; ///< Suppress outflow of particles?
+
   Field3D wall_potential; ///< Voltage of the wall. Normalised units.
 };
 
diff --git a/src/sheath_boundary_simple.cxx b/src/sheath_boundary_simple.cxx
index fde1ce56..83d88881 100644
--- a/src/sheath_boundary_simple.cxx
+++ b/src/sheath_boundary_simple.cxx
@@ -96,6 +96,11 @@ SheathBoundarySimple::SheathBoundarySimple(std::string name, Options& alloptions
           .doc("Always set phi field? Default is to only modify if already set")
           .withDefault<bool>(false);
 
+  suppress_outflow =
+      options["suppress_outflow"]
+          .doc("Suppress outflow of particles?")
+          .withDefault<bool>(false);
+
   const Options& units = alloptions["units"];
   const BoutReal Tnorm = units["eV"];
 
@@ -516,8 +521,13 @@ void SheathBoundarySimple::transform(Options& state) {
           }
 
           // Set boundary conditions on flows
-          Vi[im] = 2. * visheath - Vi[i];
-          NVi[im] = 2. * Mi * nisheath * visheath - NVi[i];
+          if (suppress_outflow == false) {
+            Vi[im] = 2. * visheath - Vi[i];
+            NVi[im] = 2. * Mi * nisheath * visheath - NVi[i];
+          } else {
+            Vi[im] = -Vi[i];
+            NVi[im] = -NVi[i];
+          }
 
           // Take into account the flow of energy due to fluid flow
           // This is additional energy flux through the sheath
@@ -580,8 +590,14 @@ void SheathBoundarySimple::transform(Options& state) {
           }
 
           // Set boundary conditions on flows
-          Vi[ip] = 2. * visheath - Vi[i];
-          NVi[ip] = 2. * Mi * nisheath * visheath - NVi[i];
+          if (suppress_outflow == false) {
+            Vi[ip] = 2. * visheath - Vi[i];
+            NVi[ip] = 2. * Mi * nisheath * visheath - NVi[i];
+          } else {
+            Vi[ip] = -Vi[i];
+            NVi[ip] = -NVi[i];
+          }
+
 
           // Take into account the flow of energy due to fluid flow
           // This is additional energy flux through the sheath
@@ -606,6 +622,11 @@ void SheathBoundarySimple::transform(Options& state) {
                       * (0.5*(coord->dx[i] + coord->dx[im]) * 0.5*(coord->dz[i] + coord->dz[im]));  // W
 
           ion_sheath_power_ylow[ip] += heatflow;       // Upper Y, so power placed in first guard cell
+
+          output << "\n****************************\n";
+          output << "visheath = " << visheath << "\n";
+          output << "ion heatflow = " << heatflow;
+          output << "\n****************************\n";
         }
       }