Merge pull request #190 from bendudson/fix-boundaries

WIP: Boundaries and misc fixes

Author: bendudson

src/electron_force_balance.cxx

@@ -18,21 +18,8 @@ void ElectronForceBalance::transform(Options &state) {
   // Note: The pressure boundary can be set in sheath boundary condition
   //       which depends on the electron velocity being set here first.
   Options& electrons = state["species"]["e"];
-  Field3D Pe = getNoBoundary<Field3D>(electrons["pressure"]);
-
-  // Note: only parallel boundaries needed
-  for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) {
-    for (int jz = 0; jz < mesh->LocalNz; jz++) {
-      Pe(r.ind, mesh->ystart - 1, jz) = Pe(r.ind, mesh->ystart, jz);
-    }
-  }
-  for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) {
-    for (int jz = 0; jz < mesh->LocalNz; jz++) {
-      Pe(r.ind, mesh->yend + 1, jz) = Pe(r.ind, mesh->yend, jz);
-    }
-  }
-
-  Field3D Ne = getNoBoundary<Field3D>(electrons["density"]);
+  Field3D Pe = GET_VALUE(Field3D, electrons["pressure"]);
+  Field3D Ne = GET_NOBOUNDARY(Field3D, electrons["density"]);
 
   ASSERT1(get<BoutReal>(electrons["charge"]) == -1.0);
 
@@ -41,8 +28,8 @@ void ElectronForceBalance::transform(Options &state) {
 
   if (IS_SET(electrons["momentum_source"])) {
     // Balance other forces from e.g. collisions
-    // Note: marked as final so can't be set later
-    force_density += get<Field3D>(electrons["momentum_source"]);
+    // Note: marked as final so can't be changed later
+    force_density += GET_VALUE(Field3D, electrons["momentum_source"]);
   }
   const Field3D Epar = force_density / floor(Ne, 1e-5);
 
@@ -54,11 +41,11 @@ void ElectronForceBalance::transform(Options &state) {
     }
     Options& species = allspecies[kv.first]; // Note: Need non-const
 
-    if (!(species.isSet("density") and species.isSet("charge"))) {
+    if (!(IS_SET(species["density"]) and IS_SET(species["charge"]))) {
       continue; // Needs both density and charge to experience a force
     }
 
-    const Field3D N = getNoBoundary<Field3D>(species["density"]);
+    const Field3D N = GET_NOBOUNDARY(Field3D, species["density"]);
     const BoutReal charge = get<BoutReal>(species["charge"]);
 
     add(species["momentum_source"],

src/evolve_energy.cxx

@@ -144,7 +144,6 @@ void EvolveEnergy::transform(Options& state) {
       P[i] = 0.0;
     }
   }
-
   P.applyBoundary("neumann");
 
   if (neumann_boundary_average_z) {
@@ -200,6 +199,37 @@ void EvolveEnergy::finally(const Options& state) {
   T = get<Field3D>(species["temperature"]);
   N = get<Field3D>(species["density"]);
   const Field3D V = get<Field3D>(species["velocity"]);
+  const BoutReal AA = get<BoutReal>(species["AA"]);
+
+  // Update boundaries of E from boundaries of P and V
+  E = toFieldAligned(E);
+  for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) {
+    for (int jz = 0; jz < mesh->LocalNz; jz++) {
+      auto i = indexAt(N, r.ind, mesh->ystart, jz);
+      auto im = i.ym();
+
+      const BoutReal psheath = 0.5 * (P[im] + P[i]);
+      const BoutReal nsheath = 0.5 * (N[im] + N[i]);
+      const BoutReal vsheath = 0.5 * (V[im] + V[i]);
+
+      const BoutReal Esheath = Cv * psheath + 0.5 * AA * nsheath * SQ(vsheath);
+      E[im] = 2 * Esheath - E[i];
+    }
+  }
+  for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) {
+    for (int jz = 0; jz < mesh->LocalNz; jz++) {
+      auto i = indexAt(N, r.ind, mesh->yend, jz);
+      auto ip = i.yp();
+
+      const BoutReal psheath = 0.5 * (P[ip] + P[i]);
+      const BoutReal nsheath = 0.5 * (N[ip] + N[i]);
+      const BoutReal vsheath = 0.5 * (V[ip] + V[i]);
+
+      const BoutReal Esheath = Cv * psheath + 0.5 * AA * nsheath * SQ(vsheath);
+      E[ip] = 2 * Esheath - E[i];
+    }
+  }
+  E = fromFieldAligned(E);
 
   Field3D Pfloor = P;
 
@@ -306,7 +336,7 @@ void EvolveEnergy::finally(const Options& state) {
   if (species.isSet("energy_source")) {
     Se += get<Field3D>(species["energy_source"]); // For diagnostic output
   }
-  if (species.isSet("energy_source")) {
+  if (species.isSet("momentum_source")) {
     Se += V * get<Field3D>(species["momentum_source"]);
   }
   ddt(E) += Se;