Merge pull request #143 from mabarnes/energy-source-term

Energy source term

docs/src/external_sources_notes.md

@@ -7,8 +7,15 @@ source term has the form of a Maxwellian with constant temperature and
 spatially-varying amplitude
 ```math
 \begin{align}
-S_i &= A_i(r,z) \exp\left( -\frac{(v_\perp^2 + v_\parallel^2)}{T_{\mathrm{source},i}} \right) \\
-S_n &= A_n(r,z) \exp\left( -\frac{(v_\zeta^2 + v_r^2 + v_z^2)}{T_{\mathrm{source},n}} \right)
+S_i &= A_i(r,z) \frac{1}{(\pi)^{3/2} (2 T_{\mathrm{source},i} / m_i)^{3/2}} \exp\left( -\frac{(v_\perp^2 + v_\parallel^2)}{T_{\mathrm{source},i}} \right) \\
+S_n &= A_n(r,z) \frac{1}{(\pi)^{3/2} (2 T_{\mathrm{source},n} / m_n)^{3/2}} \exp\left( -\frac{(v_\zeta^2 + v_r^2 + v_z^2)}{T_{\mathrm{source},n}} \right)
+\end{align}
+```
+or in 1V simulations that do not include $v_\perp$, $v_\zeta$, $v_r$ dimensions
+```math
+\begin{align}
+S_i &= A_i(r,z) \frac{1}{sqrt{\pi} \sqrt{2 T_{\mathrm{source},i} / m_i}} \exp\left( -\frac{v_\perp^2}{T_{\mathrm{source},i}} \right) \\
+S_n &= A_n(r,z) \frac{1}{sqrt{\pi} \sqrt{2 T_{\mathrm{source},n} / m_n}} \exp\left( -\frac{v_z^2}{T_{\mathrm{source},n}} \right)
 \end{align}
 ```
 
@@ -17,7 +24,7 @@ The sources are controlled by options in the `[ion_source]` and
 setting `active = true`. The constant temperature is set with the `source_T`
 option (default is 1 for ions and $T_\mathrm{wall}$ for neutrals). The
 amplitude can be set or controlled in various ways depending on the
-`controller_type` setting, as explained in the following subsection.
+`source_type` setting, as explained in the following subsection.
 
 Note that all the settings mentioned below have values given in normalised
 units (in the same way as the settings for initial profiles, etc.).
@@ -27,8 +34,8 @@ Amplitude
 
 ### Fixed amplitude (default)
 
-When `controller_type = ""` (the default), the amplitude of the source is fixed
-in time and controled by the profile options. The profile has the form
+When `source_type = "Maxwellian"` (the default), the amplitude of the source is
+fixed in time and controled by the profile options. The profile has the form
 ```math
 A(r,z) = A_0 R(r) Z(z)
 ```
@@ -52,7 +59,7 @@ available options for either are the same, so letting `x` stand for either of
 
 ### Midpoint density controller
 
-When `controller_type = "density_midpoint"` a PI controller
+When `source_type = "density_midpoint_control"` a PI controller
 ([Wikipedia](https://en.wikipedia.org/wiki/Proportional%E2%80%93integral%E2%80%93derivative_controller))
 is used to control the ion/neutral density. The 'midpoint' for the purposes of
 this controller is the point on the grid where $r=0$ and $z=0$ (there must be a
@@ -81,7 +88,7 @@ density.
 
 ### Density profile controller
 
-When `controller_type = "density_profile"` a PI controller
+When `source_type = "density_profile_control"` a PI controller
 ([Wikipedia](https://en.wikipedia.org/wiki/Proportional%E2%80%93integral%E2%80%93derivative_controller))
 is used to control the ion/neutral density profile.
 
@@ -110,7 +117,7 @@ density.
 
 The source of neutrals can be set so that some fraction of the flux of ions to
 the walls is recycled into the volume of the domain as neutrals by using the
-`controller_type = "recycling"` option.
+`source_type = "recycling"` option.
 
 The profile is set up whose spatial integral is 1
 ```math
@@ -131,6 +138,31 @@ targets.
     [`moment_kinetics.external_sources.external_neutral_source_controller!`](@ref))
     to be used in 2D simulations.
 
+### Energy source
+
+When `source_type = "energy"`, rather than just adding particles with
+temperature $T_\mathrm{source},s$, the existing plasma or neutrals in the
+domain are swapped with plasma/neutrals from a Maxwellian with
+$T_\mathrm{source},s$, so that the density is unchanged, but energy is added
+(or potentially removed if the plasma/neutrals are hotter than
+$T_\mathrm{source},s$).
+```math
+\begin{align}
+S_i &= A_i(r,z) \left[ \frac{1}{(\pi)^{3/2} (2 T_{\mathrm{source},i} / m_i)^{3/2}} \exp\left( -\frac{(v_\perp^2 + v_\parallel^2)}{T_{\mathrm{source},i}} - f_i(v_\perp, v_\parallel) \right) \right] \\
+S_n &= A_n(r,z) \left[ \frac{1}{(\pi)^{3/2} (2 T_{\mathrm{source},n} / m_n)^{3/2}} \exp\left( -\frac{(v_\zeta^2 + v_r^2 + v_z^2)}{T_{\mathrm{source},n}} - f_n(v_\zeta, v_r, v_z) \right) \right]
+\end{align}
+```
+or in 1V simulations
+```math
+\begin{align}
+S_i &= A_i(r,z) \left[ \frac{1}{sqrt{\pi} \sqrt{2 T_{\mathrm{source},i} / m_i}} \exp\left( -\frac{v_\perp^2}{T_{\mathrm{source},i}} - f_i(v_\parallel) \right) \right] \\
+S_n &= A_n(r,z) \left[ \frac{1}{sqrt{\pi} \sqrt{2 T_{\mathrm{source},n} / m_n}} \exp\left( -\frac{v_z^2}{T_{\mathrm{source},n}} - f_n(v_z) \right) \right]
+\end{align}
+```
+Note that this source does not give a fixed power input (although that might be
+a nice feature to have), it just swaps plasma/neutral particles at a constant
+rate.
+
 API
 ---
 

src/continuity.jl

@@ -31,7 +31,7 @@ function continuity_equation!(dens_out, fvec_in, moments, composition, dt, spect
     end
 
     if ion_source_settings.active
-        source_amplitude = moments.charged.external_source_amplitude
+        source_amplitude = moments.charged.external_source_density_amplitude
         @loop_s_r_z is ir iz begin
             dens_out[iz,ir,is] +=
                 dt * source_amplitude[iz,ir]
@@ -72,7 +72,7 @@ function neutral_continuity_equation!(dens_out, fvec_in, moments, composition, d
     end
 
     if neutral_source_settings.active
-        source_amplitude = moments.neutral.external_source_amplitude
+        source_amplitude = moments.neutral.external_source_density_amplitude
         @loop_s_r_z is ir iz begin
             dens_out[iz,ir,is] +=
                 dt * source_amplitude[iz,ir]

src/energy_equation.jl

@@ -23,11 +23,9 @@ function energy_equation!(ppar, fvec, moments, collisions, dt, spectral, composi
     end
 
     if ion_source_settings.active
-        source_amplitude = moments.charged.external_source_amplitude
-        source_T = ion_source_settings.source_T
+        source_amplitude = moments.charged.external_source_pressure_amplitude
         @loop_s_r_z is ir iz begin
-            ppar[iz,ir,is] += dt * source_amplitude[iz,ir] *
-                              (0.5*source_T + fvec.upar[iz,ir,is]^2)
+            ppar[iz,ir,is] += dt * source_amplitude[iz,ir]
         end
     end
 
@@ -77,11 +75,9 @@ function neutral_energy_equation!(pz, fvec, moments, collisions, dt, spectral,
     end
 
     if neutral_source_settings.active
-        source_amplitude = moments.neutral.external_source_amplitude
-        source_T = neutral_source_settings.source_T
+        source_amplitude = moments.neutral.external_source_pressure_amplitude
         @loop_s_r_z isn ir iz begin
-            pz[iz,ir,isn] += dt * source_amplitude[iz,ir] *
-                             (0.5*source_T + fvec.uz_neutral[iz,ir,isn]^2)
+            pz[iz,ir,isn] += dt * source_amplitude[iz,ir]
         end
     end