Nonorthogonal with minimal differential operators MMS test.

.pip_install_for_tests.sh

@@ -19,9 +19,8 @@ do
     fi
 done
 
-# Install Hermes-3
-git clone https://github.com/boutproject/xhermes.git
-cd xhermes
-pip3 install --user .
-cd ..
+# Install xBOUT to make sure we get latest master version
+pip3 install  --user git+https://github.com/boutproject/xBOUT.git@v0.3.7
+# Install xHermes for Hermes-3 python
+pip3 install  --user git+https://github.com/boutproject/xhermes.git
 

CMakeLists.txt

@@ -357,6 +357,28 @@ if(HERMES_TESTS)
     # Have hermes_unit_tests look for dynamic libs in ./ before searching elsewhere
     set_target_properties(hermes_unit_tests PROPERTIES BUILD_RPATH "$ORIGIN")
   endif()
+
+# Method of manufactured solutions (MMS) tests of operators
+  option(HERMES_DIVOPS_MMS_TESTS "Build the differential operator MMS tests" ON)
+  if(HERMES_DIVOPS_MMS_TESTS)
+    #TEST_SOURCES=
+    add_executable(hermes_mms_operator_tests
+                   tests/mms_operator/main.cxx)
+                   #${TEST_SOURCES})
+    target_link_libraries(hermes_mms_operator_tests PRIVATE gtest bout++::bout++ hermes-3-lib)
+    add_test(NAME differential_operators_mms_orthogonal
+           WORKING_DIRECTORY tests/mms_operator/
+           COMMAND python3 orthogonal_test.py)
+    set_tests_properties(differential_operators_mms_orthogonal PROPERTIES LABELS integration)
+    add_test(NAME differential_operators_mms_nonorthogonal
+           WORKING_DIRECTORY tests/mms_operator/
+           COMMAND python3 nonorthogonal_test.py)
+    set_tests_properties(differential_operators_mms_nonorthogonal PROPERTIES LABELS integration)
+    add_test(NAME rhs_mms_neutral_mixed
+           WORKING_DIRECTORY tests/mms_operator/
+           COMMAND python3 neutral_mixed_ddt_test.py)
+    set_tests_properties(rhs_mms_neutral_mixed PROPERTIES LABELS integration)
+  endif()
 endif()
 
 # Compile-time options

include/classical_diffusion.hxx

@@ -15,6 +15,11 @@ private:
   Field3D Dn; ///< Particle diffusion coefficient
   BoutReal custom_D; ///< User-set particle diffusion coefficient override
 
+  // Flow diagnostics
+  Field3D cls_pf_perp_xlow, cls_pf_perp_ylow;
+  Field3D cls_mf_perp_xlow, cls_mf_perp_ylow;
+  Field3D cls_nef_perp_xlow, cls_nef_perp_ylow;
+  Field3D cls_tef_perp_xlow, cls_tef_perp_ylow;
   void transform_impl(GuardedOptions& state) override;
 };
 

include/div_ops.hxx

@@ -76,6 +76,16 @@ const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f,
 const Field3D Div_par_K_Grad_par_mod(const Field3D& k, const Field3D& f, Field3D& flow_ylow,
                                      bool bndry_flux = true);
 
+/*!
+ * Div ( a Grad_perp(f) ) -- ∇⊥ ( a ⋅ ∇⊥ f) -- Vorticity
+ *
+ * This version includes corrections for non-orthogonal meshes
+ * in which the g12 and g13 components can be non-zero
+ * i.e. X-Y, X-Z and Y-Z coordinates can all be non-orthogonal.
+ */
+Field3D Div_a_Grad_perp_nonorthog(const Field3D& a, const Field3D& x, Field3D& flux_xlow,
+                                  Field3D& flux_ylow);
+
 namespace FV {
 
 /// Superbee limiter

include/neutral_mixed.hxx

@@ -18,7 +18,7 @@ struct NeutralMixed : public Component {
   /// @param options  Top-level options. Settings will be taken from options[name]
   /// @param solver   Time-integration solver to be used
   NeutralMixed(const std::string& name, Options& options, Solver *solver);
-  
+
   /// Use the final simulation state to update internal state
   /// (e.g. time derivatives)
   void finally(const Options &state) override;
@@ -30,7 +30,7 @@ struct NeutralMixed : public Component {
   void precon(const Options &state, BoutReal gamma) override;
 private:
   std::string name;  ///< Species name
-  
+
   Field3D Nn, Pn, NVn; // Density, pressure and parallel momentum
   Field3D Vn; ///< Neutral parallel velocity
   Field3D Tn; ///< Neutral temperature
@@ -45,6 +45,7 @@ private:
   Field3D DnnNn, DnnPn, DnnTn, DnnNVn; ///< Used for operators
   BoutReal flux_limit; ///< Diffusive flux limit
   BoutReal diffusion_limit;    ///< Maximum diffusion coefficient
+  BoutReal neutral_lmax;
 
   bool sheath_ydown, sheath_yup;
 
@@ -53,20 +54,23 @@ private:
   BoutReal pressure_floor; ///< Minimum Pn used when dividing Pn by Nn to get Tn.
   bool freeze_low_density; ///< Freeze evolution in low density regions?
 
-
+  BoutReal collisionality_override;     ///< Rnn input for testing
+  BoutReal density_norm, pressure_norm; ///< Normalisations
+  BoutReal momentum_norm;               ///< Normalisations
 
   bool neutral_viscosity; ///< include viscosity?
   bool neutral_conduction; ///< Include heat conduction?
   bool evolve_momentum; ///< Evolve parallel momentum?
-
+  bool normalise_sources; ///< Normalise input sources?
+
   Field3D kappa_n, eta_n; ///< Neutral conduction and viscosity
 
   bool precondition {true}; ///< Enable preconditioner?
   bool precon_laplacexy {false}; ///< Use LaplaceXY?
   bool lax_flux; ///< Use Lax flux for advection terms
   std::unique_ptr<Laplacian> inv; ///< Laplacian inversion used for preconditioning
 
-  Field3D density_source, pressure_source; ///< External input source
+  Field3D density_source, pressure_source, momentum_source; ///< External input source
   Field3D Sn, Sp, Snv; ///< Particle, pressure and momentum source
   Field3D sound_speed; ///< Sound speed for use with Lax flux
 

src/classical_diffusion.cxx

@@ -1,5 +1,6 @@
 #include "classical_diffusion.hxx"
 
+#include "../include/div_ops.hxx"
 #include <bout/fv_ops.hxx>
 
 ClassicalDiffusion::ClassicalDiffusion(std::string name, Options& alloptions, Solver*)
@@ -36,29 +37,30 @@ void ClassicalDiffusion::transform_impl(GuardedOptions& state) {
   // The only term here comes from the resistive drift
 
   Field3D Ptotal = 0.0;
-  for (auto& kv : allspecies.getChildren()) {
-    const auto species = kv.second;
+  if (!(custom_D > 0)) {
+    for (auto& kv : allspecies.getChildren()) {
+      const auto species = kv.second;
 
-    if (!(species.isSet("charge") and species.isSet("pressure"))) {
-      continue; // Skip, go to next species
-    }
-    auto q = get<BoutReal>(species["charge"]);
-    if (fabs(q) < 1e-5) {
-      continue;
+      if (!(species.isSet("charge") and species.isSet("pressure"))) {
+        continue; // Skip, go to next species
+      }
+      auto q = get<BoutReal>(species["charge"]);
+      if (fabs(q) < 1e-5) {
+        continue;
+      }
+      Ptotal += GET_VALUE(Field3D, species["pressure"]);
     }
-    Ptotal += GET_VALUE(Field3D, species["pressure"]);
   }
 
-  auto electrons = allspecies["e"];
-  const auto me = get<BoutReal>(electrons["AA"]);
-  const Field3D Ne = GET_VALUE(Field3D, electrons["density"]);
-
   // Particle diffusion coefficient. Applied to all charged species
   // so that net transport is ambipolar
 
   if (custom_D > 0) {    // User-set
     Dn = custom_D;   
   } else {                  // Calculated from collisions
+    auto electrons = allspecies["e"];
+    const auto me = get<BoutReal>(electrons["AA"]);
+    const Field3D Ne = GET_VALUE(Field3D, electrons["density"]);
     const Field3D nu_e = floor(GET_VALUE(Field3D, electrons["collision_frequency"]), 1e-10);
     Dn = floor(Ptotal, 1e-5) * me * nu_e / (floor(Ne, 1e-5) * Bsq);
   }
@@ -70,40 +72,49 @@ void ClassicalDiffusion::transform_impl(GuardedOptions& state) {
 
   for (auto kv : allspecies.getChildren()) {
     GuardedOptions species = allspecies[kv.first]; // Note: Need non-const
-
-    if (!(species.isSet("charge") and species.isSet("density"))) {
-      continue; // Skip, go to next species
-    }
-    auto q = get<BoutReal>(species["charge"]);
-    if (fabs(q) < 1e-5) {
-      continue;
+    if (!(custom_D > 0)) {
+      if (!(species.isSet("charge") and species.isSet("density"))) {
+        continue; // Skip, go to next species
+      }
+      auto q = get<BoutReal>(species["charge"]);
+      if (fabs(q) < 1e-5) {
+        continue;
+      }
     }
-
     const auto N = GET_VALUE(Field3D, species["density"]);
 
-    add(species["density_source"], FV::Div_a_Grad_perp(Dn, N));
+    // add(species["density_source"], FV::Div_a_Grad_perp(Dn, N));
+    add(species["density_source"],
+        Div_a_Grad_perp_nonorthog(Dn, N, cls_pf_perp_xlow, cls_pf_perp_ylow));
 
     if (IS_SET(species["velocity"])) {
       const auto V = GET_VALUE(Field3D, species["velocity"]);
       const auto AA = GET_VALUE(BoutReal, species["AA"]);
 
-      add(species["momentum_source"], FV::Div_a_Grad_perp(Dn * AA * V, N));
+      // add(species["momentum_source"], FV::Div_a_Grad_perp(Dn * AA * V, N));
+      add(species["momentum_source"],
+          Div_a_Grad_perp_nonorthog(Dn * AA * V, N, cls_mf_perp_xlow, cls_mf_perp_ylow));
     }
 
     if (IS_SET(species["temperature"])) {
       const auto T = GET_VALUE(Field3D, species["temperature"]);
-      add(species["energy_source"], FV::Div_a_Grad_perp(Dn * (3. / 2) * T, N));
-
-      // Cross-field heat conduction
-      // kappa_perp = 2 * n * nu_ii * rho_i^2
-
-      const auto P = GET_VALUE(Field3D, species["pressure"]);
-      const auto AA = GET_VALUE(BoutReal, species["AA"]);
+      // add(species["energy_source"], FV::Div_a_Grad_perp(Dn * (3. / 2) * T, N));
+      add(species["energy_source"],
+          Div_a_Grad_perp_nonorthog(Dn * (3. / 2) * T, N, cls_nef_perp_xlow,
+                                    cls_nef_perp_ylow));
 
       // TODO: Figure out what to do with the below
       if(custom_D < 0) {
+        // Cross-field heat conduction
+        // kappa_perp = 2 * n * nu_ii * rho_i^2
+        const auto P = GET_VALUE(Field3D, species["pressure"]);
+        const auto AA = GET_VALUE(BoutReal, species["AA"]);
         const Field3D nu = floor(GET_VALUE(Field3D, species["collision_frequency"]), 1e-10);
-        add(species["energy_source"], FV::Div_a_Grad_perp(2. * floor(P, 1e-5) * nu * AA / Bsq, T));
+        // add(species["energy_source"], FV::Div_a_Grad_perp(2. * floor(P, 1e-5) * nu * AA
+        // / Bsq, T));
+        add(species["energy_source"],
+            Div_a_Grad_perp_nonorthog(2. * floor(P, 1e-5) * nu * AA / Bsq, T,
+                                      cls_tef_perp_xlow, cls_tef_perp_ylow));
       }
     }
   }