Merge branch 'master' into neutral-advection

CMakeLists.txt

@@ -43,7 +43,8 @@ message(STATUS "Git revision: ${HERMES_REVISION}")
 option(HERMES_BUILD_BOUT "Build BOUT++ in external/BOUT-dev" ON)
 if(HERMES_BUILD_BOUT)
   set(BOUT_BUILD_EXAMPLES, OFF) # Don't create example makefiles
-  add_subdirectory(external/BOUT-dev)
+  set(HERMES_BOUT_SRC "external/BOUT-dev" CACHE STRING "Directory of BOUT++ dir")
+  add_subdirectory(${HERMES_BOUT_SRC} ${CMAKE_CURRENT_BINARY_DIR}/external/BOUT-dev)
 else()
   find_package(bout++ REQUIRED)
 endif()
@@ -224,6 +225,7 @@ if(HERMES_TESTS)
   hermes_add_integrated_test(sod-shock)
   hermes_add_integrated_test(sod-shock-energy)
   hermes_add_integrated_test(drift-wave)
+  hermes_add_integrated_test(alfven-wave)
 
   # Unit tests
   option(HERMES_UNIT_TESTS "Build the unit tests" ON)
@@ -264,5 +266,8 @@ endif()
 configure_file(include/hermes_build_config.hxx.in include/hermes_build_config.hxx)
 
 # Installation
-
-install(TARGETS hermes-3)
+if (BUILD_SHARED_LIBS)
+  install(TARGETS hermes-3 hermes-3-lib)
+else()
+  install(TARGETS hermes-3)
+endif()

docs/sphinx/code_structure.rst

@@ -54,7 +54,7 @@ for new variables which are added to the state:
     * `charge`  Charge, in units of proton charge (i.e. electron=-1)
 
     * `density`
-    * `momentum`
+    * `momentum` Parallel momentum
     * `pressure`
     * `velocity` Parallel velocity
     * `temperature`
@@ -74,12 +74,14 @@ for new variables which are added to the state:
 * `fields`
 
   * `vorticity`
-  * `phi`       Electrostatic potential
-  * `DivJdia`   Divergence of diamagnetic current
-  * `DivJcol`   Divergence of collisional current
-  * `DivJextra` Divergence of current, including 2D parallel current
-    closures.  Not including diamagnetic, parallel current due to
-    flows, or polarisation currents
+  * `phi`           Electrostatic potential
+  * `Apar`          Electromagnetic potential b dot A in induction terms
+  * `Apar_flutter`  The electromagnetic potential (b dot A) in flutter terms
+  * `DivJdia`       Divergence of diamagnetic current
+  * `DivJcol`       Divergence of collisional current
+  * `DivJextra`     Divergence of current, including 2D parallel current
+                    closures.  Not including diamagnetic, parallel current due to
+                    flows, or polarisation currents
 
 For example to get the electron density::
 

docs/sphinx/components.rst

@@ -969,7 +969,8 @@ listed after all the species groups in the component list, so that all
 the species are present in the state.
 
 One of the most important is the `collisions`_ component. This sets collision
-times for all species, which are then used 
+times for all species, which are then used in other components to calculate
+quantities like heat diffusivities and viscosity closures.
 
 .. _sound_speed:
 
@@ -1176,12 +1177,18 @@ species :math:`b` due to temperature differences, is given by:
 
 - Ion-neutral and electron-neutral collisions
 
+  *Note*: These are disabled by default. If enabled, care is needed to
+  avoid double-counting collisions in atomic reactions e.g charge-exchange
+  reactions.
+
   The cross-section for elastic collisions between charged and neutral
   particles can vary significantly. Here for simplicity we just take
   a value of :math:`5\times 10^{-19}m^2` from the NRL formulary.
 
 - Neutral-neutral collisions
 
+  *Note* This is enabled by default.
+
   The cross-section is given by
 
 .. math::
@@ -1916,6 +1923,20 @@ This functionality is not yet currently implemented for helium or neon reactions
 | R_multiplier          | Impurity species | Fixed frac. impurity radiation rate   |
 +-----------------------+------------------+---------------------------------------+
 
+The charge exchange reaction can also be modified so that the momentum transfer channel is disabled. This can be useful when
+testing the impact of the full neutral momentum equation equation compared to purely diffusive neutrals. A diffusive only model 
+leads to all of the ion momentum being lost during charge exchange due to the lack of a neutral momentum equation.
+Enabling neutral momentum introduces a more accurate transport model but also prevents CX momentum from being lost, which
+can have a significant impact on the solution and may be difficult to analyse.
+Disabling the momentum transfer channel allows you to study the impact of the improved transport only and is set as:
+
+.. code-block:: ini
+
+   [hermes]
+   components = ..., c, ...
+
+   [reactions]
+   no_neutral_cx_mom_gain = true
 
 Electromagnetic fields
 ----------------------
@@ -2017,9 +2038,26 @@ the potential in time as a diffusion equation.
 .. doxygenstruct:: RelaxPotential
    :members:
 
+.. _electromagnetic:
+
 electromagnetic
 ~~~~~~~~~~~~~~~
 
+**Notes**: When using this module,
+
+1. Set ``sound_speed:alfven_wave=true`` so that the shear Alfven wave
+   speed is included in the calculation of the fastest parallel wave
+   speed for numerical dissipation.
+2. For tokamak simulations use Neumann boundary condition on the core
+   and Dirichlet on SOL and PF boundaries by setting
+   ``electromagnetic:apar_core_neumann=true`` (this is the default).
+3. Set the potential core boundary to be constant in Y by setting
+   ``vorticity:phi_core_averagey = true``
+4. Magnetic flutter terms must be enabled to be active
+   (``electromagnetic:magnetic_flutter=true``).  They use an
+   ``Apar_flutter`` field, not the ``Apar`` field that is calculated
+   from the induction terms.
+
 This component modifies the definition of momentum of all species, to
 include the contribution from the electromagnetic potential
 :math:`A_{||}`.
@@ -2032,8 +2070,17 @@ Assumes that "momentum" :math:`p_s` calculated for all species
    p_s = m_s n_s v_{||s} + Z_s e n_s A_{||}
 
 which arises once the electromagnetic contribution to the force on
-each species is included in the momentum equation. This is normalised
-so that in dimensionless quantities
+each species is included in the momentum equation. This requires
+an additional term in the momentum equation:
+
+.. math::
+
+   \frac{\partial p_s}{\partial t} = \cdots + Z_s e A_{||} \frac{\partial n_s}{\partial t}
+
+This is implemented so that the density time-derivative is calculated using the lowest order
+terms (parallel flow, ExB drift and a low density numerical diffusion term).
+
+The above equations are normalised so that in dimensionless quantities:
 
 .. math::
 
@@ -2068,5 +2115,22 @@ To convert the species momenta into a current, we take the sum of
 
    - \frac{1}{\beta_{em}} \nabla_\perp^2 A_{||} + \sum_s \frac{Z^2 n_s}{A}A_{||} = \sum_s \frac{Z}{A} p_s
 
+The toroidal variation of density :math:`n_s` must be kept in this
+equation.  By default the iterative "Naulin" solver is used to do
+this: A fast FFT-based method is used in a fixed point iteration,
+correcting for the density variation.
+
+Magnetic flutter terms are disabled by default, and can be enabled by setting
+
+.. code-block:: ini
+
+   [electromagnetic]
+   magnetic_flutter = true
+
+This writes an ``Apar_flutter`` field to the state, which then enables perturbed
+parallel derivative terms in the ``evolve_density``, ``evolve_pressure``, ``evolve_energy`` and
+``evolve_momentum`` components. Parallel flow terms are modified, and parallel heat
+conduction.
+
 .. doxygenstruct:: Electromagnetic
    :members:

docs/sphinx/figs/alfven-wave.png

@@ -0,0 +1 @@
+../../../tests/integrated/alfven-wave/alfven-wave.png

docs/sphinx/figs/drift-wave.png

@@ -0,0 +1 @@
+../../../tests/integrated/drift-wave/drift-wave.png

docs/sphinx/tests.rst

@@ -203,3 +203,54 @@ where
 This is a cubic dispersion relation, so we find the three roots (using
 NumPy), and choose the root with the most positive growth rate
 (imaginary component of :math:`\omega`).
+
+.. figure:: figs/drift-wave.png
+   :name: drift-wave
+   :alt: Comparison of drift-wave growth rate (top) and frequency (bottom)
+   :width: 60%
+
+Alfven wave
+-----------
+
+The equations solved are
+
+.. math::
+
+   \begin{aligned}
+   \frac{\partial}{\partial t}\nabla\cdot\left(\frac{n_0 m_i}{B^2}\nabla_\perp\phi\right) =& \nabla_{||}J_{||} = -\nabla_{||}\left(en_ev_{||e}\right) \\
+   \frac{\partial}{\partial t}\left(m_en_ev_{||e} - en_eA_{||}\right) =& -\nabla\cdot\left(m_en_ev_{||e} \mathbf{b}v_{||e}\right) + en_e\partial_{||}\phi - 0.51\nu_{ei}n_im_ev_{||e} \\
+   J_{||} =& \frac{1}{\mu_0}\nabla_\perp^2 A_{||}
+   \end{aligned}
+
+Linearising around a stationary background with constant density
+:math:`n_0` and temperature :math:`T_0`, using
+:math:`\frac{\partial}{\partial t}\rightarrow -i\omega` gives:
+
+.. math::
+
+   \begin{aligned}
+   \tilde{\phi} =& -\frac{k_{||}}{\omega k_\perp^2}\frac{eB^2}{m_i}\tilde{v_{||e}} \\
+   \omega \left( m_e \tilde{v_{||e}} - e\tilde{A}_{||}\right) =& -ek_{||}\tilde{\phi} - i0.51\nu_{ei}m_e\tilde{v_{||e}} \\
+   en_0\tilde{v_{||e}} =& -\frac{k_\perp^2}{\mu_0}\tilde{A}_{||}
+   \end{aligned}
+
+Rearranging results in a quadratic dispersion relation:
+
+.. math::
+
+   \omega^2\left(1 + \frac{k_\perp^2 c^2}{\omega_{pe}^2}\right) + i 0.51\nu_{ei}\frac{k_\perp^2 c^2}{\omega_{pe}^2}\omega - k_{||}^2V_A^2 = 0
+
+where :math:`V_A = B / \sqrt{\mu_0 n_0 m_i}` is the Alfven speed, and
+:math:`c / \omega_{pe} = \sqrt{m_e / \left(\mu_0 n_0 e^2\right)}` is
+the electron skin depth.
+
+When collisions are neglected, we obtain the result
+
+.. math::
+
+   \omega^2 = V_A^2\frac{k_{||}^2}{1 + k_\perp^2 c^2 / \omega_{pe}^2}
+
+.. figure:: figs/alfven-wave.png
+   :name: alfven-wave
+   :alt: Alfven wave speed, as function of parallel and perpendicular wavenumbers
+   :width: 60%

include/amjuel_reaction.hxx

@@ -164,7 +164,7 @@ protected:
         Ne.getRegion("RGN_NOBNDRY"))(Ne, N1, Te);
 
     // Loss is reduced by heating
-    energy_loss -= (electron_heating / Tnorm) * reaction_rate;
+    energy_loss -= (electron_heating / Tnorm) * reaction_rate * radiation_multiplier;
 
     subtract(electron["energy_source"], energy_loss);
   }

include/binormal_stpm.hxx

@@ -31,12 +31,15 @@ struct BinormalSTPM : public Component {
   ///     - density correction
   ///
   void transform(Options& state) override;
+  void outputVars(Options &state) override;
 
 
 private:
   std::string name; ///< Short name of the species e.g. h+
+  bool diagnose; ///< Output diagnostics?
   Field3D Theta, chi, D, nu; ///< Field line pitch, anomalous thermal, momentum diffusion
   Field3D nu_Theta, chi_Theta, D_Theta; ///< nu/Theta, chi/Theta, D/Theta, precalculated
+  Field3D Theta_inv; ///< Precalculate 1/Theta
 
 };
 

include/component.hxx

@@ -319,4 +319,15 @@ void set_with_attrs(Options& option, T value, std::initializer_list<std::pair<st
   option.setAttributes(attrs);
 }
 
+#if CHECKLEVEL >= 1
+template<>
+inline void set_with_attrs(Options& option, Field3D value, std::initializer_list<std::pair<std::string, Options::AttributeType>> attrs) {
+  if (!value.isAllocated()) {
+    throw BoutException("set_with_attrs: Field3D assigned to {:s} is not allocated", option.str());
+  }
+  option.force(value);
+  option.setAttributes(attrs);
+}
+#endif
+
 #endif // HERMES_COMPONENT_H

include/div_ops.hxx

@@ -23,8 +23,8 @@
 
  */
 
-#ifndef HERMES_DIV_OPS_H
-#define HERMES_DIV_OPS_H
+#ifndef DIV_OPS_H
+#define DIV_OPS_H
 
 #include <bout/field3d.hxx>
 #include <bout/fv_ops.hxx>
@@ -44,6 +44,11 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f,
                                     bool bndry_flux = true, bool poloidal = false,
                                     bool positive = false);
 
+/// This version has an extra coefficient 'g' that is linearly interpolated
+/// onto cell faces
+const Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D &n, const Field3D &g, const Field3D &f, 
+                                    bool bndry_flux = true, bool positive = false);
+
 const Field3D Div_Perp_Lap_FV_Index(const Field3D& a, const Field3D& f, bool xflux);
 
 const Field3D Div_Z_FV_Index(const Field3D& a, const Field3D& f);
@@ -210,9 +215,10 @@ const Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in,
             flux = bndryval * vpar * vpar;
           } else {
             // Add flux due to difference in boundary values
-            flux = s.R * vpar * sv.R
-                   + BOUTMAX(wave_speed(i, j, k), fabs(v(i, j, k)), fabs(v(i, j + 1, k)))
-                         * (s.R * sv.R - bndryval * vpar);
+            flux =
+                s.R * sv.R * sv.R
+                + BOUTMAX(wave_speed(i, j, k), fabs(v(i, j, k)), fabs(v(i, j + 1, k)))
+                  * (s.R * sv.R - bndryval * vpar);
           }
         } else {
           // Maximum wave speed in the two cells
@@ -238,9 +244,10 @@ const Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in,
             flux = bndryval * vpar * vpar;
           } else {
             // Add flux due to difference in boundary values
-            flux = s.L * vpar * sv.L
-                   - BOUTMAX(wave_speed(i, j, k), fabs(v(i, j, k)), fabs(v(i, j - 1, k)))
-                         * (s.L * sv.L - bndryval * vpar);
+            flux =
+                s.L * sv.L * sv.L
+                - BOUTMAX(wave_speed(i, j, k), fabs(v(i, j, k)), fabs(v(i, j - 1, k)))
+                  * (s.L * sv.L - bndryval * vpar);
           }
         } else {
           // Maximum wave speed in the two cells
@@ -683,4 +690,4 @@ const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Fi
 
 } // namespace FV
 
-#endif //  HERMES_DIV_OPS_H
+#endif //  DIV_OPS_H

include/electromagnetic.hxx

@@ -64,6 +64,13 @@ private:
 
   std::unique_ptr<Laplacian> aparSolver; // Laplacian solver in X-Z
 
+  bool const_gradient; // Set neumann boundaries by extrapolation
+  BoutReal apar_boundary_timescale; // Relaxation timescale
+  BoutReal last_time;  // The last time the boundaries were updated
+
+  bool magnetic_flutter; ///< Set the magnetic flutter term?
+  Field3D Apar_flutter;
+
   bool diagnose; ///< Output additional diagnostics?
 };
 

include/evolve_momentum.hxx

@@ -35,7 +35,8 @@ private:
   std::string name;     ///< Short name of species e.g "e"
 
   Field3D NV;           ///< Species parallel momentum (normalised, evolving)
-  Field3D NV_solver;    ///< Momentum as input from solver
+  Field3D NV_err;       ///< Difference from momentum as input from solver
+  Field3D NV_solver;    ///< Momentum as calculated in the solver
   Field3D V;            ///< Species parallel velocity
 
   Field3D momentum_source; ///< From other components. Stored for diagnostic output

include/fixed_fraction_radiation.hxx

@@ -302,7 +302,7 @@ struct FixedFractionRadiation : public Component {
     AUTO_TRACE();
 
     if (diagnose) {
-      set_with_attrs(state[std::string("R") + name], radiation,
+      set_with_attrs(state[std::string("R") + name], -radiation,
                      {{"time_dimension", "t"},
                       {"units", "W / m^3"},
                       {"conversion", SI::qe * Tnorm * Nnorm * FreqNorm},