Merge pull request #95 from bendudson/paper-inputs

WIP: Inputs for Hermes-3 paper

docs/sphinx/components.rst

@@ -650,10 +650,52 @@ term:
 
 .. math::
 
-   Q_{ab,F} = - F_{ab} u_a
+   Q_{ab,F} = \frac{m_b}{m_a + m_b} \left( u_b - u_a \right) F_{ab}
+
+This term has some important properties:
+
+1. It is always positive: Collisions of two species with the same
+   temperature never leads to cooling.
+2. It is Galilean invariant: Shifting both species' velocity by the
+   same amount leaves :math:`Q_{ab,F}` unchanged.
+3. If both species have the same mass, the thermal energy
+   change due to slowing down is shared equally between them.
+4. If one species is much heavier than the other, for example
+   electron-ion collisions, the lighter species is preferentially
+   heated. This recovers e.g. Braginskii expressions for :math:`Q_{ei}`
+   and :math:`Q_{ie}`.
+
+This can be derived by considering the exchange of energy
+:math:`W_{ab,F}` between two species at the same temperature but
+different velocities. If the pressure is evolved then it contains
+a term that balances the change in kinetic energy due to changes
+in velocity:
 
-Energy exchange, heat transferred to species `a` from species `b` due to temperature
-differences, is given by:
+.. math::
+
+   \begin{aligned}
+   \frac{\partial}{\partial t}\left(m_a n_a u_a\right) =& \ldots + F_{ab} \\
+   \frac{\partial}{\partial t}\left(\frac{3}{2}p_a\right) =& \ldots - F_{ab} u_a + W_{ab, F}
+   \end{aligned}
+
+For momentum and energy conservation we must have :math:`F_{ab}=-F_{ba}`
+and :math:`W_{ab,F} = -W_{ba,F}`. Comparing the above to the
+`Braginskii expression
+<https://farside.ph.utexas.edu/teaching/plasma/lectures/node35.html>`_
+we see that for ion-electron collisions the term :math:`- F_{ab}u_a + W_{ab, F}`
+goes to zero, so :math:`W_{ab, F} \sim u_aF_{ab}` for
+:math:`m_a \gg m_b`. An expression that has all these desired properties
+is
+
+.. math::
+
+   W_{ab,F} = \left(\frac{m_a u_a + m_b u_a}{m_a + m_b}\right)F_{ab}
+
+which is not Galilean invariant but when combined with the :math:`- F_{ab} u_a`
+term gives a change in pressure that is invariant, as required.
+
+Thermal energy exchange, heat transferred to species :math:`a` from
+species :math:`b` due to temperature differences, is given by:
 
 .. math::
 
@@ -817,6 +859,56 @@ Notes:
    The reason for this convention is the existence of the inverse reactions:
    `t + d+ -> t+ + d` outputs diagnostics `Ftd+_cx` and `Fd+t_cx`.
 
+2. Reactions typically convert species from one to another, leading to
+   a transfer of mass momentum and energy. For a reaction converting
+   species :math:`a` to species :math:`b` at rate :math:`R` (units
+   of events per second per volume) we have transfers:
+
+   .. math::
+
+      \begin{aligned}
+      \frac{\partial}{\partial t} n_a =& \ldots - R \\
+      \frac{\partial}{\partial t} n_b =& \ldots + R \\
+      \frac{\partial}{\partial t}\left( m n_a u_a\right) =& \ldots + F_{ab} \\
+      \frac{\partial}{\partial t}\left( m n_a u_a\right) =& \ldots + F_{ba} \\
+      \frac{\partial}{\partial t}\left( \frac{3}{2} p_a \right) =& \ldots - F_{ab}u_a + W_{ab} - \frac{1}{2}mRu_a^2 \\
+      \frac{\partial}{\partial t}\left( \frac{3}{2} p_b \right) =& \ldots - F_{ba}u_b + W_{ba} + \frac{1}{2}mRu_b^2
+      \end{aligned}
+      
+  where both species have the same mass: :math:`m_a = m_b = m`. In the
+  pressure equations the :math:`-F_{ab}u_a` comes from splitting the
+  kinetic and thermal energies; :math:`W_{ab}=-W_{ba}` is the energy
+  transfer term that we need to find; The final term balances the loss
+  of kinetic energy at fixed momentum due to a particle source or
+  sink.
+
+  The momentum transfer :math:`F_{ab}=-F{ba}` is the momentum carried
+  by the converted ions: :math:`F_{ab}=-m R u_a`. To find
+  :math:`W_{ab}` we note that for :math:`p_a = 0` the change in pressure
+  must go to zero: :math:`-F_{ab}u_a + W_{ab} -\frac{1}{2}mRu_a^2 = 0`.
+
+  .. math::
+
+      \begin{aligned}
+      W_{ab} =& F_{ab}u_a + \frac{1}{2}mRu_a^2 \\
+      =& - mR u_a^2 + \frac{1}{2}mRu_a^2\\
+      =& -\frac{1}{2}mRu_a^2
+      \end{aligned}
+
+  Substituting into the above gives:
+
+  .. math::
+
+     \begin{aligned}
+     \frac{\partial}{\partial t}\left( \frac{3}{2} p_b \right) =& \ldots - F_{ba}u_b + W_{ba} + \frac{1}{2}mRu_b^2 \\
+     =& \ldots - mRu_au_b + \frac{1}{2}mRu_a^2 + \frac{1}{2}mRu_a^2 \\
+     =& \ldots + \frac{1}{2}mR\left(u_a - u_b\right)^2
+     \end{aligned}
+
+  This has the property that the change in pressure of both species is
+  Galilean invariant. This transfer term is included in the Amjuel reactions
+  and hydrogen charge exchange.
+
 Hydrogen
 ~~~~~~~~
 

examples/1D-recycling/BOUT.inp

@@ -13,14 +13,14 @@
 #  Does not include recombination
 
 
-nout = 200
-timestep = 2000
+nout = 400
+timestep = 5000
 
 MXG = 0  # No guard cells in X
 
 [mesh]
 nx = 1
-ny = 200   # Resolution along field-line
+ny = 400   # Resolution along field-line
 nz = 1
 
 length = 30           # Length of the domain in meters
@@ -53,7 +53,13 @@ Bnorm = 1
 Tnorm = 100
 
 [solver]
-mxstep = 100000
+type = beuler  # Backward Euler steady-state solver
+snes_type = newtonls
+ksp_type = gmres
+max_nonlinear_iterations = 10
+
+atol = 1e-7
+rtol = 1e-5
 
 [sheath_boundary]
 
@@ -77,13 +83,16 @@ charge = 1
 AA = 2
 
 density_upstream = 1e19  # Upstream density [m^-3]
+density_source_positive = false  # Force source to be > 0?
+density_controller_i = 5e-4
+density_controller_p = 1e-2
 
 thermal_conduction = true  # in evolve_pressure
 
 diagnose = true
 
 recycle_as = d
-recycle_multiplier = 0.99  # Recycling fraction
+recycle_multiplier = 1  # Recycling fraction
 
 [Nd+]
 
@@ -148,6 +157,7 @@ species = d+
 
 [reactions]
 type = (
-        d + e -> d+ + 2e,   # Deuterium ionisation
-        d + d+ -> d+ + d,   # Charge exchange
+        d + e -> d+ + 2e,     # Deuterium ionisation
+        d+ + e -> d,          # Deuterium recombination
+        d + d+ -> d+ + d,     # Charge exchange
        )

examples/1D-recycling/makeplot.py

@@ -46,6 +46,7 @@
 
     plt.title(f"Time {i}")
     plt.savefig(f"1d_recycling_{i:02d}.png")
+    plt.savefig(f"1d_recycling_{i:02d}.pdf")
     plt.close()
 
 os.system('convert -resize 50% -delay 20 -loop 0 *.png 1d_recycling.gif')

examples/1D-recycling/plot_convergence.py

@@ -0,0 +1,65 @@
+# Analyse convergence towards steady state
+
+paths = [".", "pvode"]
+labels = ["NK", "PVODE"]
+linestyles = ["-", "--"]
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from boutdata import collect
+
+cumulative_calls = []
+times = []
+source_mags = []
+
+for path in paths:
+    try:
+        # Number of calls
+        ncalls = collect("ncalls", path=path)  # number of calls each output
+
+        cumulative_calls.append(np.cumsum(ncalls))   # Summed over outputs
+
+        t_array = collect("t_array", path=path)
+        wci = collect("Omega_ci", path=path)
+        times.append(t_array / wci) # Seconds
+
+        # density source
+        m = collect("density_source_multiplier_d+", path=path)
+        source_mags.append(np.abs(m + 1))
+    except:
+        break
+
+for calls, source, label, linestyle in zip(cumulative_calls, source_mags, labels, linestyles):
+    plt.plot(calls, source, label=label, linestyle=linestyle)
+
+plt.yscale('log')
+plt.xlabel('RHS evaluations')
+plt.ylabel('Source magnitude [arb]')
+plt.legend()
+plt.savefig("source_rhsevals.pdf")
+plt.show()
+
+for time, source, label, linestyle in zip(times, source_mags, labels, linestyles):
+    plt.plot(time, source, label=label, linestyle=linestyle)
+plt.yscale('log')
+plt.xlabel('Time [s]')
+plt.ylabel('Source magnitude [arb]')
+plt.legend()
+plt.savefig("source_time.pdf")
+plt.show()
+
+# Time derivatives
+
+for path, calls, label, linestyle in zip(paths, cumulative_calls, labels, linestyles):
+    for varname, color in [("ddt(Nd+)", 'k'), ("ddt(Pd+)", 'r'), ("ddt(NVd+)",'b')]:
+        dv = collect(varname, path=path)
+        dv_rms = np.sqrt(np.sum(dv[:,0,:,0]**2, axis=1))
+        plt.plot(calls, dv_rms, label=label + " " + varname, linestyle=linestyle, color=color)
+plt.yscale('log')
+plt.ylabel("RMS time derivative over domain")
+plt.xlabel("RHS evaluations")
+plt.legend()
+plt.savefig("timederivs_rhsevals.pdf")
+plt.show()
+