SecondaryElectronEnergyBC

Kinetic secondary electron for mean electron energy boundary condition

Overview

SecondaryElectronEnergyBC is an electron mean energy density of secondary electrons induced by ion flux outflow boundary condition. SecondaryElectronEnergyBC assumes the electrostatic approximation for the electric field.

$var element = document.getElementById("moose-equation-3289f19a-7536-45ba-9e5e-d6e9feacf2aa");katex.render("a_{e} = \\begin{cases} 1, & \\mu_{e} \\left( \\nabla (V) \\right) \\cdot \\textbf{n} > 0\\\\ 0, & \\mu_{e} \\left( \\nabla (V) \\right)\\cdot \\textbf{n} \\leq 0\\\\ \\end{cases} \\\\[10pt] a_{i} = \\begin{cases} 1, & \\mu_{i} \\left( - \\nabla (V) \\right) \\cdot \\textbf{n} > 0\\\\ 0, & \\mu_{i} \\left( - \\nabla (V) \\right)\\cdot \\textbf{n} \\leq 0\\\\ \\end{cases} \\\\[10pt] \\Gamma_{\\epsilon} \\cdot \\textbf{n} = \\frac{1 - r_{e}}{1 + r_{e}} \\frac{\\text{-}5}{6} v_{\\text{th},e} n_{\\gamma} \\varepsilon_{\\gamma} - \\frac{5}{3} \\frac{2}{1+r} (1 - a_{e}) \\gamma \\Gamma_{j} \\cdot \\textbf{n} \\\\[10pt] \\Gamma_{j} \\cdot \\textbf{n} = \\frac{1-r_{j}}{1+r_{j}} \\left[ (2 a_{j}-1) \\ \\mu_{j} \\left( \\text{-} \\nabla (V) \\right) \\cdot \\textbf{n} \\ n_{j} + \\frac{1}{2}v_{\\text{th},j}n_{j} \\right] \\\\[10pt] v_{\\text{th},e} = \\sqrt{\\frac{8e}{\\pi m_{e}} \\frac{2}{3} \\frac{n_{\\varepsilon}}{n_{e}}} \\\\[10pt] v_{\\text{th},j} = \\sqrt{\\frac{8k_{B}T}{\\pi m_{j}}} \\\\[10pt] n_{\\gamma} = (1 - a_{e}) \\frac{\\gamma \\Gamma_{i} \\cdot \\textbf{n}}{\\mu_{e} (\\text{-}\\nabla V) \\cdot \\textbf{n}}", element, {displayMode:true,throwOnError:false});$

Where the subscripts $var element = document.getElementById("moose-equation-aa895f5d-7e93-48ae-8b4b-1a77b637b6f5");katex.render("e", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-a5e4f9b1-06f6-4c9c-ba41-dd558a559e3f");katex.render("j", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-7c920370-d7cb-4280-880b-c705f7a93d48");katex.render("\\varepsilon", element, {displayMode:false,throwOnError:false});$ represents properties of electrons, ions and electron energy respectively, $var element = document.getElementById("moose-equation-16c5afbb-a86f-464c-9c47-bd96da9ffdde");katex.render("\\Gamma", element, {displayMode:false,throwOnError:false});$ is the flux of the electron or ions, $var element = document.getElementById("moose-equation-1e306b0c-3f2d-438d-ac8e-1b1e7041afa4");katex.render("\\textbf{n}", element, {displayMode:false,throwOnError:false});$ is the normal vector of the boundary, $var element = document.getElementById("moose-equation-331b1b26-3e68-49ad-b8d8-3d0be4bb88a3");katex.render("\\mu", element, {displayMode:false,throwOnError:false});$ is the mobility coefficient, $var element = document.getElementById("moose-equation-fbe6d1cd-e721-4c0a-9a4b-1cf4e7558ddf");katex.render("n", element, {displayMode:false,throwOnError:false});$ is the species density, $var element = document.getElementById("moose-equation-6fe279a4-3231-4cb7-a8d7-47ce25302627");katex.render("v_\\text{th}", element, {displayMode:false,throwOnError:false});$ is the thermal velocity of the species, $var element = document.getElementById("moose-equation-fc2f1c38-094c-445b-bbb7-392e62338c0f");katex.render("k_{B}", element, {displayMode:false,throwOnError:false});$ is the Boltzmann constant, $var element = document.getElementById("moose-equation-042134a3-fa43-476a-93c5-19e8ca238904");katex.render("T", element, {displayMode:false,throwOnError:false});$ is the gas temperature, $var element = document.getElementById("moose-equation-45a26c95-51da-462a-ac11-2aa507fb031a");katex.render("n_{\\gamma}", element, {displayMode:false,throwOnError:false});$ is electron density emitted by the surface, and $var element = document.getElementById("moose-equation-2ef308b6-ce3c-47db-99a0-363f108c1a08");katex.render("\\varepsilon_{\\gamma}", element, {displayMode:false,throwOnError:false});$ is the energy of the secondary electrons. $var element = document.getElementById("moose-equation-aa0dbe62-b6c5-4a60-867d-4cd03b4311e9");katex.render("a", element, {displayMode:false,throwOnError:false});$ is defined such that the outflow is non-zero when the drift velocity is directed towards the wall and zero otherwise. $var element = document.getElementById("moose-equation-aa8218c1-261b-4d02-bfa2-19a4e8ff9dc7");katex.render("r", element, {displayMode:false,throwOnError:false});$ is defined as the fraction of particles reflected by the surface. When converting the density to log form and applying a scaling factor of the mesh, the strong form for SecondaryElectronEnergyBC is defined as

$var element = document.getElementById("moose-equation-b525d804-63b2-4d80-8929-1f839f1525e0");katex.render("v_{\\text{th},e} = \\sqrt{\\frac{8e}{\\pi m_{e}} \\frac{2}{3} \\exp (N_{\\varepsilon} - N_{e})} \\\\[10pt] \\Gamma_{j} \\cdot \\textbf{n} = \\frac{1-r_{j}}{1+r_{j}} \\left[ (2 a_{j}-1) \\ \\mu_{j} \\left( \\text{-} \\nabla (V / l_{c}) \\right) \\cdot \\textbf{n} \\ \\exp(N_{j}) + \\frac{1}{2}v_{\\text{th},j}\\exp(N_{j}) \\right]", element, {displayMode:true,throwOnError:false});$

Where $var element = document.getElementById("moose-equation-f8edfca0-ed39-4bbf-8794-3c60846f5c89");katex.render("N_{e}", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-5557432f-ccbd-4fe7-9ea8-6ad30462a80c");katex.render("N_{i}", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-b690b8d2-3aae-45c5-8f32-5dcbc90982bb");katex.render("N_{\\varepsilon}", element, {displayMode:false,throwOnError:false});$ is the molar density of the electrons, ions and electron energy in logarithmic form and $var element = document.getElementById("moose-equation-a2462209-905c-420f-9a5b-3186488a6ac0");katex.render("l_{c}", element, {displayMode:false,throwOnError:false});$ is the scaling factor of the mesh.

Example Input File Syntax

[BCs]
  [secondary_energy_left]
    type = SecondaryElectronEnergyBC
    variable = mean_en
    boundary = 'left'
    potential = potential
    electrons = em
    ions = 'Arp'
    r = 0
    emission_coeffs = 0.05
    secondary_electron_energy = 3
    position_units = ${dom0Scale}
  []
[]

Input Parameters

boundaryThe list of boundary IDs from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Unit:(no unit assumed)
Controllable:No
Description:The list of boundary IDs from the mesh where this object applies
electronsThe electron density in log form
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The electron density in log form
emission_coeffsA species-dependent list of secondary electron emission coefficients
C++ Type:std::vector<std::string>
Unit:(no unit assumed)
Controllable:No
Description:A species-dependent list of secondary electron emission coefficients
ionsA list of ion densities in log form
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:A list of ion densities in log form
position_unitsUnits of position.
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Units of position.
potentialThe electric potential
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The electric potential
rThe reflection coefficient of the electrons.
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The reflection coefficient of the electrons.
secondary_electron_energyThe secondary electron energy in eV
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The secondary electron energy in eV
variableThe name of the variable that this residual object operates on
C++ Type:NonlinearVariableName
Unit:(no unit assumed)
Controllable:No
Description:The name of the variable that this residual object operates on

Required Parameters

displacementsThe displacements
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The displacements
prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
r_ion0The reflection coefficient of the ions.
Default:0
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The reflection coefficient of the ions.
use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.

Optional Parameters

absolute_value_vector_tagsThe tags for the vectors this residual object should fill with the absolute value of the residual contribution
C++ Type:std::vector<TagName>
Unit:(no unit assumed)
Controllable:No
Description:The tags for the vectors this residual object should fill with the absolute value of the residual contribution
extra_matrix_tagsThe extra tags for the matrices this Kernel should fill
C++ Type:std::vector<TagName>
Unit:(no unit assumed)
Controllable:No
Description:The extra tags for the matrices this Kernel should fill
extra_vector_tagsThe extra tags for the vectors this Kernel should fill
C++ Type:std::vector<TagName>
Unit:(no unit assumed)
Controllable:No
Description:The extra tags for the vectors this Kernel should fill
matrix_tagssystemThe tag for the matrices this Kernel should fill
Default:system
C++ Type:MultiMooseEnum
Unit:(no unit assumed)
Options:nontime, system
Controllable:No
Description:The tag for the matrices this Kernel should fill
vector_tagsnontimeThe tag for the vectors this Kernel should fill
Default:nontime
C++ Type:MultiMooseEnum
Unit:(no unit assumed)
Options:nontime, time
Controllable:No
Description:The tag for the vectors this Kernel should fill

Tagging Parameters

control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Unit:(no unit assumed)
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
diag_save_inThe name of auxiliary variables to save this BC's diagonal jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Unit:(no unit assumed)
Controllable:No
Description:The name of auxiliary variables to save this BC's diagonal jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Unit:(no unit assumed)
Controllable:Yes
Description:Set the enabled status of the MooseObject.
implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
save_inThe name of auxiliary variables to save this BC's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Unit:(no unit assumed)
Controllable:No
Description:The name of auxiliary variables to save this BC's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Unit:(no unit assumed)
Controllable:No
Description:The seed for the master random number generator
use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.

Advanced Parameters

Input Files

(test/tests/automatic_differentiation/ad_argon.i)
(test/tests/DriftDiffusionAction/mean_en_no_actions.i)
(tutorial/tutorial05-PlasmaWaterInterface/DC_argon-With-Water.i)
(test/tests/1d_dc/mean_en.i)
(test/tests/DriftDiffusionAction/mean_en_actions.i)
(test/tests/1d_dc/mean_en_multi.i)
(test/tests/crane_action/townsend_units.i)

(test/tests/1d_dc/mean_en.i)

dom0Scale = 1e-3
dom1Scale = 1e-7

[GlobalParams]
  offset = 20
  # offset = 0
  potential_units = kV
  use_moles = true
  # potential_units = V
[]

[Mesh]
  [file]
    type = FileMeshGenerator
    file = 'liquidNew.msh'
  []
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = '0'
    paired_block = '1'
    new_boundary = 'master0_interface'
    input = file
  []
  [interface_again]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = '1'
    paired_block = '0'
    new_boundary = 'master1_interface'
    input = interface
  []
  [left]
    type = SideSetsFromNormalsGenerator
    normals = '-1 0 0'
    new_boundary = 'left'
    input = interface_again
  []
  [right]
    type = SideSetsFromNormalsGenerator
    normals = '1 0 0'
    new_boundary = 'right'
    input = left
  []
[]

[Problem]
  type = FEProblem
  # kernel_coverage_check = false
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
    ksp_norm = none
  []
[]

[Executioner]
  type = Transient
  end_time = 1e-1
  petsc_options = '-snes_converged_reason -snes_linesearch_monitor'
  #petsc_options = '-snes_converged_reason -snes_linesearch_monitor -snes_test_jacobian'
  # petsc_options = '-snes_test_display'
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount -ksp_type -snes_linesearch_minlambda'
  petsc_options_value = 'lu NONZERO 1.e-10 preonly 1e-3'
  # petsc_options_iname = '-pc_type -sub_pc_type'
  # petsc_options_value = 'asm lu'
  # petsc_options_iname = '-snes_type'
  # petsc_options_value = 'test'
  nl_rel_tol = 1e-4
  nl_abs_tol = 7.6e-5
  dtmin = 1e-12
  l_max_its = 20
  [TimeSteppers]
    [Adaptive]
      type = IterationAdaptiveDT
      cutback_factor = 0.4
      dt = 1e-11
      # dt = 1.1
      growth_factor = 1.2
      optimal_iterations = 15
    []
  []
[]

[Outputs]
  perf_graph = true
  # print_linear_residuals = false
  [out]
    type = Exodus
    execute_on = 'final'
  []
  [dof_map]
    type = DOFMap
  []
[]

[Debug]
  show_var_residual_norms = true
[]

[UserObjects]
  [data_provider]
    type = ProvideMobility
    electrode_area = 5.02e-7 # Formerly 3.14e-6
    ballast_resist = 1e6
    e = 1.6e-19
  []
[]

[Kernels]
  [em_time_deriv]
    type = ElectronTimeDerivative
    variable = em
    block = 0
  []
  [em_advection]
    type = EFieldAdvection
    variable = em
    potential = potential
    block = 0
    position_units = ${dom0Scale}
  []
  [em_diffusion]
    type = CoeffDiffusion
    variable = em
    block = 0
    position_units = ${dom0Scale}
  []
  [em_ionization]
    type = ElectronsFromIonization
    variable = em
    potential = potential
    mean_en = mean_en
    em = em
    block = 0
    position_units = ${dom0Scale}
  []
  [em_log_stabilization]
    type = LogStabilizationMoles
    variable = em
    block = 0
  []

  [emliq_time_deriv]
    type = ElectronTimeDerivative
    variable = emliq
    block = 1
  []
  [emliq_advection]
    type = EFieldAdvection
    variable = emliq
    potential = potential
    block = 1
    position_units = ${dom1Scale}
  []
  [emliq_diffusion]
    type = CoeffDiffusion
    variable = emliq
    block = 1
    position_units = ${dom1Scale}
  []
  [emliq_reactant_first_order_rxn]
    type = ReactantFirstOrderRxn
    variable = emliq
    block = 1
  []
  [emliq_water_bi_sink]
    type = ReactantAARxn
    variable = emliq
    block = 1
  []
  [emliq_log_stabilization]
    type = LogStabilizationMoles
    variable = emliq
    block = 1
  []

  [potential_diffusion_dom1]
    type = CoeffDiffusionLin
    variable = potential
    block = 0
    position_units = ${dom0Scale}
  []
  [potential_diffusion_dom2]
    type = CoeffDiffusionLin
    variable = potential
    block = 1
    position_units = ${dom1Scale}
  []
  [Arp_charge_source]
    type = ChargeSourceMoles_KV
    variable = potential
    charged = Arp
    block = 0
  []
  [em_charge_source]
    type = ChargeSourceMoles_KV
    variable = potential
    charged = em
    block = 0
  []
  [emliq_charge_source]
    type = ChargeSourceMoles_KV
    variable = potential
    charged = emliq
    block = 1
  []
  [OHm_charge_source]
    type = ChargeSourceMoles_KV
    variable = potential
    charged = OHm
    block = 1
  []

  [Arp_time_deriv]
    type = ElectronTimeDerivative
    variable = Arp
    block = 0
  []
  [Arp_advection]
    type = EFieldAdvection
    variable = Arp
    potential = potential
    position_units = ${dom0Scale}
    block = 0
  []
  [Arp_diffusion]
    type = CoeffDiffusion
    variable = Arp
    block = 0
    position_units = ${dom0Scale}
  []
  [Arp_ionization]
    type = IonsFromIonization
    variable = Arp
    potential = potential
    em = em
    mean_en = mean_en
    block = 0
    position_units = ${dom0Scale}
  []
  [Arp_log_stabilization]
    type = LogStabilizationMoles
    variable = Arp
    block = 0
  []

  [OHm_time_deriv]
    type = ElectronTimeDerivative
    variable = OHm
    block = 1
  []
  [OHm_advection]
    type = EFieldAdvection
    variable = OHm
    potential = potential
    block = 1
    position_units = ${dom1Scale}
  []
  [OHm_diffusion]
    type = CoeffDiffusion
    variable = OHm
    block = 1
    position_units = ${dom1Scale}
  []
  [OHm_log_stabilization]
    type = LogStabilizationMoles
    variable = OHm
    block = 1
  []
  [OHm_product_first_order_rxn]
    type = ProductFirstOrderRxn
    variable = OHm
    v = emliq
    block = 1
  []
  [OHm_product_aabb_rxn]
    type = ProductAABBRxn
    variable = OHm
    v = emliq
    block = 1
  []

  [mean_en_time_deriv]
    type = ElectronTimeDerivative
    variable = mean_en
    block = 0
  []
  [mean_en_advection]
    type = EFieldAdvection
    variable = mean_en
    potential = potential
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_diffusion]
    type = CoeffDiffusion
    variable = mean_en
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_joule_heating]
    type = JouleHeating
    variable = mean_en
    potential = potential
    em = em
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_ionization]
    type = ElectronEnergyLossFromIonization
    variable = mean_en
    potential = potential
    em = em
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_elastic]
    type = ElectronEnergyLossFromElastic
    variable = mean_en
    potential = potential
    em = em
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_excitation]
    type = ElectronEnergyLossFromExcitation
    variable = mean_en
    potential = potential
    em = em
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_log_stabilization]
    type = LogStabilizationMoles
    variable = mean_en
    block = 0
    offset = 15
  []
[]

[Variables]
  [potential]
  []
  [em]
    block = 0
  []
  [emliq]
    block = 1
    # scaling = 1e-5
  []

  [Arp]
    block = 0
  []

  [mean_en]
    block = 0
    # scaling = 1e-1
  []

  [OHm]
    block = 1
    # scaling = 1e-5
  []
[]

[AuxVariables]
  [e_temp]
    block = 0
    order = CONSTANT
    family = MONOMIAL
  []
  [x]
    order = CONSTANT
    family = MONOMIAL
  []
  [x_node]
  []
  [rho]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [rholiq]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [em_lin]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [emliq_lin]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [Arp_lin]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [OHm_lin]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [Efield]
    order = CONSTANT
    family = MONOMIAL
  []
  [Current_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [Current_emliq]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [Current_Arp]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [Current_OHm]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [tot_gas_current]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [tot_liq_current]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [tot_flux_OHm]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [EFieldAdvAux_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [DiffusiveFlux_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [EFieldAdvAux_emliq]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [DiffusiveFlux_emliq]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [PowerDep_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [PowerDep_Arp]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [ProcRate_el]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [ProcRate_ex]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [ProcRate_iz]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
[]

[AuxKernels]
  [PowerDep_em]
    type = ADPowerDep
    density_log = em
    potential = potential
    art_diff = false
    potential_units = kV
    variable = PowerDep_em
    position_units = ${dom0Scale}
    block = 0
  []
  [PowerDep_Arp]
    type = ADPowerDep
    density_log = Arp
    potential = potential
    art_diff = false
    potential_units = kV
    variable = PowerDep_Arp
    position_units = ${dom0Scale}
    block = 0
  []
  [ProcRate_el]
    type = ADProcRate
    em = em
    potential = potential
    proc = el
    variable = ProcRate_el
    position_units = ${dom0Scale}
    block = 0
  []
  [ProcRate_ex]
    type = ADProcRate
    em = em
    potential = potential
    proc = ex
    variable = ProcRate_ex
    position_units = ${dom0Scale}
    block = 0
  []
  [ProcRate_iz]
    type = ADProcRate
    em = em
    potential = potential
    proc = iz
    variable = ProcRate_iz
    position_units = ${dom0Scale}
    block = 0
  []
  [e_temp]
    type = ElectronTemperature
    variable = e_temp
    electron_density = em
    mean_en = mean_en
    block = 0
  []
  [x_g]
    type = Position
    variable = x
    position_units = ${dom0Scale}
    block = 0
  []
  [x_l]
    type = Position
    variable = x
    position_units = ${dom1Scale}
    block = 1
  []
  [x_ng]
    type = Position
    variable = x_node
    position_units = ${dom0Scale}
    block = 0
  []
  [x_nl]
    type = Position
    variable = x_node
    position_units = ${dom1Scale}
    block = 1
  []
  [rho]
    type = ParsedAux
    variable = rho
    coupled_variables = 'em_lin Arp_lin'
    expression = 'Arp_lin - em_lin'
    execute_on = 'timestep_end'
    block = 0
  []
  [rholiq]
    type = ParsedAux
    variable = rholiq
    coupled_variables = 'emliq_lin OHm_lin' # H3Op_lin OHm_lin'
    expression = '-emliq_lin - OHm_lin' # 'H3Op_lin - em_lin - OHm_lin'
    execute_on = 'timestep_end'
    block = 1
  []
  [tot_gas_current]
    type = ParsedAux
    variable = tot_gas_current
    coupled_variables = 'Current_em Current_Arp'
    expression = 'Current_em + Current_Arp'
    execute_on = 'timestep_end'
    block = 0
  []
  [tot_liq_current]
    type = ParsedAux
    variable = tot_liq_current
    coupled_variables = 'Current_emliq Current_OHm' # Current_H3Op Current_OHm'
    expression = 'Current_emliq + Current_OHm' # + Current_H3Op + Current_OHm'
    execute_on = 'timestep_end'
    block = 1
  []
  [em_lin]
    type = DensityMoles
    variable = em_lin
    density_log = em
    block = 0
  []
  [emliq_lin]
    type = DensityMoles
    variable = emliq_lin
    density_log = emliq
    block = 1
  []
  [Arp_lin]
    type = DensityMoles
    variable = Arp_lin
    density_log = Arp
    block = 0
  []
  [OHm_lin]
    type = DensityMoles
    variable = OHm_lin
    density_log = OHm
    block = 1
  []
  [Efield_g]
    type = Efield
    component = 0
    potential = potential
    variable = Efield
    position_units = ${dom0Scale}
    block = 0
  []
  [Efield_l]
    type = Efield
    component = 0
    potential = potential
    variable = Efield
    position_units = ${dom1Scale}
    block = 1
  []
  [Current_em]
    type = ADCurrent
    potential = potential
    density_log = em
    variable = Current_em
    art_diff = false
    block = 0
    position_units = ${dom0Scale}
  []
  [Current_emliq]
    type = ADCurrent
    potential = potential
    density_log = emliq
    variable = Current_emliq
    art_diff = false
    block = 1
    position_units = ${dom1Scale}
  []
  [Current_Arp]
    type = ADCurrent
    potential = potential
    density_log = Arp
    variable = Current_Arp
    art_diff = false
    block = 0
    position_units = ${dom0Scale}
  []
  [Current_OHm]
    block = 1
    type = ADCurrent
    potential = potential
    density_log = OHm
    variable = Current_OHm
    art_diff = false
    position_units = ${dom1Scale}
  []
  [tot_flux_OHm]
    block = 1
    type = ADTotalFlux
    potential = potential
    density_log = OHm
    variable = tot_flux_OHm
  []
  [EFieldAdvAux_em]
    type = ADEFieldAdvAux
    potential = potential
    density_log = em
    variable = EFieldAdvAux_em
    block = 0
    position_units = ${dom0Scale}
  []
  [DiffusiveFlux_em]
    type = ADDiffusiveFlux
    density_log = em
    variable = DiffusiveFlux_em
    block = 0
    position_units = ${dom0Scale}
  []
  [EFieldAdvAux_emliq]
    type = ADEFieldAdvAux
    potential = potential
    density_log = emliq
    variable = EFieldAdvAux_emliq
    block = 1
    position_units = ${dom1Scale}
  []
  [DiffusiveFlux_emliq]
    type = ADDiffusiveFlux
    density_log = emliq
    variable = DiffusiveFlux_emliq
    block = 1
    position_units = ${dom1Scale}
  []
[]

[InterfaceKernels]
  [em_advection]
    type = InterfaceAdvection
    mean_en_neighbor = mean_en
    potential_neighbor = potential
    neighbor_var = em
    variable = emliq
    boundary = master1_interface
    position_units = ${dom1Scale}
    neighbor_position_units = ${dom0Scale}
  []
  [em_diffusion]
    type = InterfaceLogDiffusionElectrons
    mean_en_neighbor = mean_en
    neighbor_var = em
    variable = emliq
    boundary = master1_interface
    position_units = ${dom1Scale}
    neighbor_position_units = ${dom0Scale}
  []
[]

[BCs]
  [mean_en_physical_right]
    type = HagelaarEnergyBC
    variable = mean_en
    boundary = 'master0_interface'
    potential = potential
    electrons = em
    r = 0.99
    position_units = ${dom0Scale}
  []
  [mean_en_physical_left]
    type = HagelaarEnergyBC
    variable = mean_en
    boundary = 'left'
    potential = potential
    electrons = em
    r = 0
    position_units = ${dom0Scale}
  []
  [secondary_energy_left]
    type = SecondaryElectronEnergyBC
    variable = mean_en
    boundary = 'left'
    potential = potential
    electrons = em
    ions = 'Arp'
    r = 0
    emission_coeffs = 0.05
    secondary_electron_energy = 3
    position_units = ${dom0Scale}
  []

  [potential_left]
    type = NeumannCircuitVoltageMoles_KV
    variable = potential
    boundary = left
    function = potential_bc_func
    ions = Arp
    data_provider = data_provider
    electrons = em
    electron_energy = mean_en
    r = 0
    position_units = ${dom0Scale}
    emission_coeffs = 0.05
  []
  [potential_dirichlet_right]
    type = DirichletBC
    variable = potential
    boundary = right
    value = 0
    preset = false
  []
  [em_physical_right]
    type = HagelaarElectronBC
    variable = em
    boundary = 'master0_interface'
    potential = potential
    electron_energy = mean_en
    r = 0.99
    position_units = ${dom0Scale}
  []
  [Arp_physical_right_diffusion]
    type = HagelaarIonDiffusionBC
    variable = Arp
    boundary = 'master0_interface'
    r = 0
    position_units = ${dom0Scale}
  []
  [Arp_physical_right_advection]
    type = HagelaarIonAdvectionBC
    variable = Arp
    boundary = 'master0_interface'
    potential = potential
    r = 0
    position_units = ${dom0Scale}
  []

  [em_physical_left]
    type = HagelaarElectronBC
    variable = em
    boundary = 'left'
    potential = potential
    electron_energy = mean_en
    r = 0
    position_units = ${dom0Scale}
  []
  [sec_electrons_left]
    type = SecondaryElectronBC
    variable = em
    boundary = 'left'
    potential = potential
    ions = Arp
    electron_energy = mean_en
    r = 0
    position_units = ${dom0Scale}
    emission_coeffs = 0.05
  []
  [Arp_physical_left_diffusion]
    type = HagelaarIonDiffusionBC
    variable = Arp
    boundary = 'left'
    r = 0
    position_units = ${dom0Scale}
  []
  [Arp_physical_left_advection]
    type = HagelaarIonAdvectionBC
    variable = Arp
    boundary = 'left'
    potential = potential
    r = 0
    position_units = ${dom0Scale}
  []

  [emliq_right]
    type = DCIonBC
    variable = emliq
    boundary = right
    potential = potential
    position_units = ${dom1Scale}
  []
  [OHm_physical]
    type = DCIonBC
    variable = OHm
    boundary = 'right'
    potential = potential
    position_units = ${dom1Scale}
  []
[]

[ICs]
  [em_ic]
    type = ConstantIC
    variable = em
    value = -21
    block = 0
  []
  [emliq_ic]
    type = ConstantIC
    variable = emliq
    value = -21
    block = 1
  []
  [Arp_ic]
    type = ConstantIC
    variable = Arp
    value = -21
    block = 0
  []
  [mean_en_ic]
    type = ConstantIC
    variable = mean_en
    value = -20
    block = 0
  []
  [OHm_ic]
    type = ConstantIC
    variable = OHm
    value = -15.6
    block = 1
  []
  [potential_ic]
    type = FunctionIC
    variable = potential
    function = potential_ic_func
  []
[]

[Functions]
  [potential_bc_func]
    type = ParsedFunction
    # expression = '1.25*tanh(1e6*t)'
    expression = -1.25
  []
  [potential_ic_func]
    type = ParsedFunction
    expression = '-1.25 * (1.0001e-3 - x)'
  []
[]

[Materials]
  [gas_block]
    type = Gas
    interp_trans_coeffs = true
    interp_elastic_coeff = true
    ramp_trans_coeffs = false
    em = em
    potential = potential
    ip = Arp
    mean_en = mean_en
    user_se_coeff = 0.05
    block = 0
    property_tables_file = td_argon_mean_en.txt
  []
  [water_block]
    type = Water
    block = 1
    potential = potential
  []
[]

(test/tests/automatic_differentiation/ad_argon.i)

dom0Scale = 1e-3

[GlobalParams]
  offset = 30
  potential_units = kV
  use_moles = true
[]

[Mesh]
  [file]
    type = FileMeshGenerator
    file = 'ad_argon.msh'
  []
  [left]
    type = SideSetsFromNormalsGenerator
    normals = '-1 0 0'
    new_boundary = 'left'
    input = file
  []
  [right]
    type = SideSetsFromNormalsGenerator
    normals = '1 0 0'
    new_boundary = 'right'
    input = left
  []
[]

[Problem]
  type = FEProblem
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  automatic_scaling = true
  compute_scaling_once = false
  end_time = 1e-1
  petsc_options = '-snes_converged_reason -snes_linesearch_monitor'
  solve_type = NEWTON
  line_search = 'basic'
  petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount'
  petsc_options_value = 'lu NONZERO 1.e-10'
  nl_rel_tol = 1e-6
  nl_abs_tol = 1e-10
  dtmin = 1e-18
  l_max_its = 20
  [TimeSteppers]
    [Adaptive]
      type = IterationAdaptiveDT
      cutback_factor = 0.4
      dt = 1e-11
      growth_factor = 1.2
      optimal_iterations = 30
    []
  []
[]

[Outputs]
  perf_graph = true
  [out]
    type = Exodus
    execute_on = 'final'
  []
[]

[UserObjects]
  [data_provider]
    type = ProvideMobility
    electrode_area = 5.02e-7 # Formerly 3.14e-6
    ballast_resist = 1e6
    e = 1.6e-19
  []
[]

[Kernels]
  [Arex_time_deriv]
    type = TimeDerivativeLog
    variable = Ar*
    block = 0
  []
  [Arex_diffusion]
    type = CoeffDiffusion
    variable = Ar*
    block = 0
    position_units = ${dom0Scale}
  []
  [Arex_log_stabilization]
    type = LogStabilizationMoles
    variable = Ar*
    block = 0
    #offset = 25
  []

  [em_time_deriv]
    type = TimeDerivativeLog
    variable = em
    block = 0
  []
  [em_advection]
    type = EFieldAdvection
    variable = em
    potential = potential
    block = 0
    position_units = ${dom0Scale}
  []
  [em_diffusion]
    type = CoeffDiffusion
    variable = em
    block = 0
    position_units = ${dom0Scale}
  []

  [em_log_stabilization]
    type = LogStabilizationMoles
    variable = em
    block = 0
  []

  [potential_diffusion_dom1]
    type = CoeffDiffusionLin
    variable = potential
    block = 0
    position_units = ${dom0Scale}
  []
  [Arp_charge_source]
    type = ChargeSourceMoles_KV
    variable = potential
    charged = Arp
    block = 0
  []
  [Ar2p_charge_source]
    type = ChargeSourceMoles_KV
    variable = potential
    charged = Ar2p
    block = 0
  []
  [em_charge_source]
    type = ChargeSourceMoles_KV
    variable = potential
    charged = em
    block = 0
  []

  [Arp_time_deriv]
    type = TimeDerivativeLog
    variable = Arp
    block = 0
  []
  [Arp_advection]
    type = EFieldAdvection
    variable = Arp
    potential = potential
    position_units = ${dom0Scale}
    block = 0
  []
  [Arp_diffusion]
    type = CoeffDiffusion
    variable = Arp
    block = 0
    position_units = ${dom0Scale}
  []

  [Arp_log_stabilization]
    type = LogStabilizationMoles
    variable = Arp
    block = 0
  []

  [Ar2p_time_deriv]
    #type = ElectronTimeDerivative
    type = TimeDerivativeLog
    variable = Ar2p
    block = 0
  []
  [Ar2p_advection]
    type = EFieldAdvection
    variable = Ar2p
    potential = potential
    position_units = ${dom0Scale}
    block = 0
  []
  [Ar2p_diffusion]
    type = CoeffDiffusion
    variable = Ar2p
    block = 0
    position_units = ${dom0Scale}
  []

  [Ar2p_log_stabilization]
    type = LogStabilizationMoles
    variable = Ar2p
    block = 0
  []

  [mean_en_time_deriv]
    type = TimeDerivativeLog
    variable = mean_en
    block = 0
  []
  [mean_en_advection]
    type = EFieldAdvection
    variable = mean_en
    potential = potential
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_diffusion]
    type = CoeffDiffusion
    variable = mean_en
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_joule_heating]
    type = JouleHeating
    variable = mean_en
    potential = potential
    em = em
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_log_stabilization]
    type = LogStabilizationMoles
    variable = mean_en
    block = 0
    offset = 15
  []
[]

[Variables]
  [potential]
  []

  [em]
    initial_condition = -24
    block = 0
  []

  [Arp]
    initial_condition = -24.69314718056
    block = 0
  []

  [Ar*]
    initial_condition = -26
    block = 0
  []

  [Ar2p]
    initial_condition = -24.69314718056
    block = 0
  []

  [mean_en]
    block = 0
    initial_condition = -24
  []
[]

[AuxVariables]
  [Arex_lin]
    block = 0
    order = CONSTANT
    family = MONOMIAL
  []
  [Ar2p_lin]
    block = 0
    order = CONSTANT
    family = MONOMIAL
  []
  [Ar]
    block = 0
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 3.70109
  []
  [e_temp]
    block = 0
    order = CONSTANT
    family = MONOMIAL
  []
  [x]
    order = CONSTANT
    family = MONOMIAL
  []
  [x_node]
  []
  [rho]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [em_lin]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [Arp_lin]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [Efield]
    order = CONSTANT
    family = MONOMIAL
  []
  [Current_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [Current_Arp]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [Current_Ar2p]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [tot_gas_current]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [EFieldAdvAux_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [DiffusiveFlux_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [PowerDep_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [PowerDep_Arp]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
[]

[AuxKernels]
  [Arex_lin]
    type = DensityMoles
    variable = Arex_lin
    density_log = Ar*
    block = 0
  []
  [Ar2p_lin]
    type = DensityMoles
    variable = Ar2p_lin
    density_log = Ar2p
    block = 0
  []
  [PowerDep_em]
    type = ADPowerDep
    density_log = em
    potential = potential
    art_diff = false
    potential_units = kV
    variable = PowerDep_em
    position_units = ${dom0Scale}
    block = 0
  []
  [PowerDep_Arp]
    type = ADPowerDep
    density_log = Arp
    potential = potential
    art_diff = false
    potential_units = kV
    variable = PowerDep_Arp
    position_units = ${dom0Scale}
    block = 0
  []
  [e_temp]
    type = ElectronTemperature
    variable = e_temp
    electron_density = em
    mean_en = mean_en
    block = 0
  []
  [x_g]
    type = Position
    variable = x
    position_units = ${dom0Scale}
    block = 0
  []
  [x_ng]
    type = Position
    variable = x_node
    position_units = ${dom0Scale}
    block = 0
  []
  [rho]
    type = ParsedAux
    variable = rho
    coupled_variables = 'em_lin Arp_lin Ar2p_lin'
    expression = 'Arp_lin + Ar2p_lin - em_lin'
    execute_on = 'timestep_end'
    block = 0
  []
  [tot_gas_current]
    type = ParsedAux
    variable = tot_gas_current
    coupled_variables = 'Current_em Current_Arp Current_Ar2p'
    expression = 'Current_em + Current_Arp + Current_Ar2p'
    execute_on = 'timestep_end'
    block = 0
  []
  [em_lin]
    type = DensityMoles
    variable = em_lin
    density_log = em
    block = 0
  []
  [Arp_lin]
    type = DensityMoles
    variable = Arp_lin
    density_log = Arp
    block = 0
  []
  [Efield_g]
    type = Efield
    component = 0
    potential = potential
    variable = Efield
    position_units = ${dom0Scale}
    block = 0
  []
  [Current_em]
    type = ADCurrent
    potential = potential
    density_log = em
    variable = Current_em
    art_diff = false
    block = 0
    position_units = ${dom0Scale}
  []
  [Current_Arp]
    type = ADCurrent
    potential = potential
    density_log = Arp
    variable = Current_Arp
    art_diff = false
    block = 0
    position_units = ${dom0Scale}
  []
  [Current_Ar2p]
    type = ADCurrent
    potential = potential
    density_log = Ar2p
    variable = Current_Ar2p
    art_diff = false
    block = 0
    position_units = ${dom0Scale}
  []
  [EFieldAdvAux_em]
    type = ADEFieldAdvAux
    potential = potential
    density_log = em
    variable = EFieldAdvAux_em
    block = 0
    position_units = ${dom0Scale}
  []
  [DiffusiveFlux_em]
    type = ADDiffusiveFlux
    density_log = em
    variable = DiffusiveFlux_em
    block = 0
    position_units = ${dom0Scale}
  []
[]

[BCs]
  [mean_en_physical_right]
    type = HagelaarEnergyBC
    variable = mean_en
    boundary = 'right'
    potential = potential
    electrons = em
    r = 0.0
    position_units = ${dom0Scale}
  []
  [mean_en_physical_left]
    type = HagelaarEnergyBC
    variable = mean_en
    boundary = 'left'
    potential = potential
    electrons = em
    r = 0
    position_units = ${dom0Scale}
  []
  [secondary_energy_left]
    type = SecondaryElectronEnergyBC
    variable = mean_en
    boundary = 'left'
    potential = potential
    electrons = em
    ions = 'Arp Ar2p'
    r = 0
    emission_coeffs = '0.05 0.05'
    position_units = ${dom0Scale}
    secondary_electron_energy = 3
  []

  [potential_left]
    type = NeumannCircuitVoltageMoles_KV
    variable = potential
    boundary = left
    function = potential_bc_func
    ions = 'Arp Ar2p'
    data_provider = data_provider
    electrons = em
    electron_energy = mean_en
    r = 0
    position_units = ${dom0Scale}
    emission_coeffs = '0.05 0.05'
  []
  [potential_dirichlet_right]
    type = DirichletBC
    variable = potential
    boundary = right
    value = 0
    preset = false
  []
  [em_physical_right]
    type = HagelaarElectronBC
    variable = em
    boundary = 'right'
    potential = potential
    electron_energy = mean_en
    r = 0.0
    position_units = ${dom0Scale}
  []
  [Arp_physical_right_diffusion]
    type = HagelaarIonDiffusionBC
    variable = Arp
    boundary = 'right'
    r = 0
    position_units = ${dom0Scale}
  []
  [Arp_physical_right_advection]
    type = HagelaarIonAdvectionBC
    variable = Arp
    boundary = 'right'
    potential = potential
    r = 0
    position_units = ${dom0Scale}
  []
  [Ar2p_physical_right_diffusion]
    type = HagelaarIonDiffusionBC
    variable = Ar2p
    boundary = 'right'
    r = 0
    position_units = ${dom0Scale}
  []
  [Ar2p_physical_right_advection]
    type = HagelaarIonAdvectionBC
    variable = Ar2p
    boundary = 'right'
    potential = potential
    r = 0
    position_units = ${dom0Scale}
  []

  [em_physical_left]
    type = HagelaarElectronBC
    variable = em
    boundary = 'left'
    potential = potential
    electron_energy = mean_en
    r = 0
    position_units = ${dom0Scale}
  []
  [sec_electrons_left]
    type = SecondaryElectronBC
    variable = em
    boundary = 'left'
    potential = potential
    ions = Arp
    electron_energy = mean_en
    r = 0
    position_units = ${dom0Scale}
    emission_coeffs = 0.05
  []
  [Arp_physical_left_diffusion]
    type = HagelaarIonDiffusionBC
    variable = Arp
    boundary = 'left'
    r = 0
    position_units = ${dom0Scale}
  []
  [Arp_physical_left_advection]
    type = HagelaarIonAdvectionBC
    variable = Arp
    boundary = 'left'
    potential = potential
    r = 0
    position_units = ${dom0Scale}
  []
  [Ar2p_physical_left_diffusion]
    type = HagelaarIonDiffusionBC
    variable = Ar2p
    boundary = 'left'
    r = 0
    position_units = ${dom0Scale}
  []
  [Ar2p_physical_left_advection]
    type = HagelaarIonAdvectionBC
    variable = Ar2p
    boundary = 'left'
    potential = potential
    r = 0
    position_units = ${dom0Scale}
  []
[]

[ICs]
  [potential_ic]
    type = FunctionIC
    variable = potential
    function = potential_ic_func
  []
[]

[Functions]
  [potential_bc_func]
    type = ParsedFunction
    expression = -0.8
  []
  [potential_ic_func]
    type = ParsedFunction
    expression = '-0.8 * (1 - x)'
  []
[]

[Materials]
  [electron_moments]
    type = ADGasElectronMoments
    block = 0
    em = em
    mean_en = mean_en
    property_tables_file = 'argon_chemistry_rates/electron_moments.txt'
  []

  [gas_constants]
    type = GenericConstantMaterial
    block = 0
    prop_names = ' e         N_A     k_boltz  eps       T_gas massem   p_gas  n_gas    se_coeff se_energy'
    prop_values = '1.6e-19 6.022e23  1.38e-23 8.854e-12 300   9.11e-31 1.01e5 40.4915  0.05     3.'
  []
  [ad_gas_constants]
    type = ADGenericConstantMaterial
    block = 0
    prop_names = 'diffpotential'
    prop_values = '8.85e-12'
  []
  [gas_species_0]
    type = ADHeavySpecies
    heavy_species_name = Arp
    heavy_species_mass = 6.64e-26
    heavy_species_charge = 1.0
    block = 0
  []
  [gas_species_1]
    type = ADHeavySpecies
    heavy_species_name = Ar2p
    heavy_species_mass = 1.328e-25
    heavy_species_charge = 1.0
    block = 0
  []

  [gas_species_2]
    type = ADHeavySpecies
    heavy_species_name = Ar
    heavy_species_mass = 6.64e-26
    heavy_species_charge = 0.0
    block = 0
  []

  [gas_species_3]
    type = ADHeavySpecies
    heavy_species_name = Ar*
    heavy_species_mass = 6.64e-26
    heavy_species_charge = 0.0
    block = 0
  []
[]

[Reactions]
  # This argon reaction network based on a ZDPlasKin example:
  # zdplaskin.laplace.univ-tlse.fr
  # Example: "Micro-cathode sustained discharged in Ar"
  [Argon]
    species = 'em Arp Ar2p Ar*'
    aux_species = 'Ar'
    reaction_coefficient_format = 'townsend'
    gas_species = 'Ar'
    electron_energy = 'mean_en'
    electron_density = 'em'
    include_electrons = true
    file_location = 'argon_chemistry_rates'
    equation_constants = 'Tgas'
    equation_values = '300'

    # For function-based reactions (e.g. those whose rate coefficients are
    # enclosed in curly braces, {}), the e_temp auxiliary variable is
    # named as a variable. Note that while e_temp is an AuxVariable, it does
    # depend on nonlinear variables and therefore should contribute to the
    # jacobian of any functions which depend on it.
    # At the moment this jacobian contribution is not included, but the
    # contribution is minimal compared to advection and diffusion so we
    # can get away with a NEWTON solver in this case.
    #
    # For more complex reaction networks with many functions that depend on
    # electron temperature, it may be necessary to use PJFNK instead of
    # NEWTON.
    equation_variables = 'e_temp'
    potential = 'potential'
    use_log = true
    position_units = ${dom0Scale}
    use_ad = true
    block = 0

    reactions = 'em + Ar -> em + Ar               : EEDF [elastic] (reaction1)
                 em + Ar -> em + Ar*              : EEDF [-11.5]   (reaction2)
                 em + Ar -> em + em + Arp         : EEDF [-15.76]  (reaction3)
                 em + Ar* -> em + Ar              : EEDF [11.5]    (reaction4)
                 em + Ar* -> em + em + Arp        : EEDF [-4.3]    (reaction5)
                 Ar2p + em -> Ar* + Ar            : {5.1187e11*(e_temp/300)^(-0.67)}
                 Ar2p + Ar -> Arp + Ar + Ar       : {3.649332e12/Tgas*exp(-15130/Tgas)}
                 Ar* + Ar* -> Ar2p + em           : 3.6132e8
                 Arp + em + em -> Ar + em         : {3.17314235e9*(e_temp/11600)^(-4.5)}
                 Ar* + Ar + Ar -> Ar + Ar + Ar    : 5077.02776
                 Arp + Ar + Ar -> Ar2p + Ar       : {81595.089*(Tgas/300)^(-0.4)}'
  []
[]

(test/tests/DriftDiffusionAction/mean_en_no_actions.i)

#This is the input file that supplied the gold output file.
#It is the same as the input file in tests/1d_dc/mean_en.i,
#execpt some of the Aux Variables are renamed for the Action test

dom0Scale = 1e-3
dom1Scale = 1e-7

[GlobalParams]
  offset = 20
  # offset = 0
  potential_units = kV
  use_moles = true
  # potential_units = V
[]

[Mesh]
  [file]
    type = FileMeshGenerator
    file = 'liquidNew.msh'
  []
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = '0'
    paired_block = '1'
    new_boundary = 'master0_interface'
    input = file
  []
  [interface_again]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = '1'
    paired_block = '0'
    new_boundary = 'master1_interface'
    input = interface
  []
  [left]
    type = SideSetsFromNormalsGenerator
    normals = '-1 0 0'
    new_boundary = 'left'
    input = interface_again
  []
  [right]
    type = SideSetsFromNormalsGenerator
    normals = '1 0 0'
    new_boundary = 'right'
    input = left
  []
[]

[Problem]
  type = FEProblem
  # kernel_coverage_check = false
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
    ksp_norm = none
  []
[]

[Executioner]
  type = Transient
  end_time = 1e-1
  petsc_options = '-snes_converged_reason -snes_linesearch_monitor'
  # petsc_options = '-snes_test_display'
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount -ksp_type -snes_linesearch_minlambda'
  petsc_options_value = 'lu NONZERO 1.e-10 preonly 1e-3'
  # petsc_options_iname = '-pc_type -sub_pc_type'
  # petsc_options_value = 'asm lu'
  # petsc_options_iname = '-snes_type'
  # petsc_options_value = 'test'
  nl_rel_tol = 1e-4
  nl_abs_tol = 7.6e-5
  dtmin = 1e-12
  l_max_its = 20
  [TimeSteppers]
    [Adaptive]
      type = IterationAdaptiveDT
      cutback_factor = 0.4
      dt = 1e-11
      # dt = 1.1
      growth_factor = 1.2
      optimal_iterations = 15
    []
  []
[]

[Outputs]
  file_base = out
  perf_graph = true
  [out]
    type = Exodus
    execute_on = 'final'
  []
[]

[Debug]
  show_var_residual_norms = true
[]

[UserObjects]
  [data_provider]
    type = ProvideMobility
    electrode_area = 5.02e-7 # Formerly 3.14e-6
    ballast_resist = 1e6
    e = 1.6e-19
    # electrode_area = 1.1
    # ballast_resist = 1.1
    # e = 1.1
  []
[]

[Kernels]
  [em_time_deriv]
    type = ElectronTimeDerivative
    variable = em
    block = 0
  []
  [em_advection]
    type = EFieldAdvection
    variable = em
    potential = potential
    block = 0
    position_units = ${dom0Scale}
  []
  [em_diffusion]
    type = CoeffDiffusion
    variable = em
    block = 0
    position_units = ${dom0Scale}
  []
  [em_ionization]
    type = ElectronsFromIonization
    variable = em
    potential = potential
    mean_en = mean_en
    em = em
    block = 0
    position_units = ${dom0Scale}
  []
  [em_log_stabilization]
    type = LogStabilizationMoles
    variable = em
    block = 0
  []
  # [em_advection_stabilization]
  #   type = EFieldArtDiff
  #   variable = em
  #   potential = potential
  #   block = 0
  # []

  [emliq_time_deriv]
    type = ElectronTimeDerivative
    variable = emliq
    block = 1
  []
  [emliq_advection]
    type = EFieldAdvection
    variable = emliq
    potential = potential
    block = 1
    position_units = ${dom1Scale}
  []
  [emliq_diffusion]
    type = CoeffDiffusion
    variable = emliq
    block = 1
    position_units = ${dom1Scale}
  []
  [emliq_reactant_first_order_rxn]
    type = ReactantFirstOrderRxn
    variable = emliq
    block = 1
  []
  [emliq_water_bi_sink]
    type = ReactantAARxn
    variable = emliq
    block = 1
  []
  [emliq_log_stabilization]
    type = LogStabilizationMoles
    variable = emliq
    block = 1
  []
  # [emliq_advection_stabilization]
  #   type = EFieldArtDiff
  #   variable = emliq
  #   potential = potential
  #   block = 1
  # []

  [potential_diffusion_dom1]
    type = CoeffDiffusionLin
    variable = potential
    block = 0
    position_units = ${dom0Scale}
  []
  [potential_diffusion_dom2]
    type = CoeffDiffusionLin
    variable = potential
    block = 1
    position_units = ${dom1Scale}
  []
  [Arp_charge_source]
    type = ChargeSourceMoles_KV
    variable = potential
    charged = Arp
    block = 0
  []
  [em_charge_source]
    type = ChargeSourceMoles_KV
    variable = potential
    charged = em
    block = 0
  []
  [emliq_charge_source]
    type = ChargeSourceMoles_KV
    variable = potential
    charged = emliq
    block = 1
  []
  [OHm_charge_source]
    type = ChargeSourceMoles_KV
    variable = potential
    charged = OHm
    block = 1
  []

  [Arp_time_deriv]
    type = ElectronTimeDerivative
    variable = Arp
    block = 0
  []
  [Arp_advection]
    type = EFieldAdvection
    variable = Arp
    potential = potential
    position_units = ${dom0Scale}
    block = 0
  []
  [Arp_diffusion]
    type = CoeffDiffusion
    variable = Arp
    block = 0
    position_units = ${dom0Scale}
  []
  [Arp_ionization]
    type = IonsFromIonization
    variable = Arp
    potential = potential
    em = em
    mean_en = mean_en
    block = 0
    position_units = ${dom0Scale}
  []
  [Arp_log_stabilization]
    type = LogStabilizationMoles
    variable = Arp
    block = 0
  []
  # [Arp_advection_stabilization]
  #   type = EFieldArtDiff
  #   variable = Arp
  #   potential = potential
  #   block = 0
  # []

  [OHm_time_deriv]
    type = ElectronTimeDerivative
    variable = OHm
    block = 1
  []
  [OHm_advection]
    type = EFieldAdvection
    variable = OHm
    potential = potential
    block = 1
    position_units = ${dom1Scale}
  []
  [OHm_diffusion]
    type = CoeffDiffusion
    variable = OHm
    block = 1
    position_units = ${dom1Scale}
  []
  [OHm_log_stabilization]
    type = LogStabilizationMoles
    variable = OHm
    block = 1
  []
  # [OHm_advection_stabilization]
  #   type = EFieldArtDiff
  #   variable = OHm
  #   potential = potential
  #   block = 1
  # []
  [OHm_product_first_order_rxn]
    type = ProductFirstOrderRxn
    variable = OHm
    v = emliq
    block = 1
  []
  [OHm_product_aabb_rxn]
    type = ProductAABBRxn
    variable = OHm
    v = emliq
    block = 1
  []

  [mean_en_time_deriv]
    type = ElectronTimeDerivative
    variable = mean_en
    block = 0
  []
  [mean_en_advection]
    type = EFieldAdvection
    variable = mean_en
    potential = potential
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_diffusion]
    type = CoeffDiffusion
    variable = mean_en
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_joule_heating]
    type = JouleHeating
    variable = mean_en
    potential = potential
    em = em
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_ionization]
    type = ElectronEnergyLossFromIonization
    variable = mean_en
    potential = potential
    em = em
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_elastic]
    type = ElectronEnergyLossFromElastic
    variable = mean_en
    potential = potential
    em = em
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_excitation]
    type = ElectronEnergyLossFromExcitation
    variable = mean_en
    potential = potential
    em = em
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_log_stabilization]
    type = LogStabilizationMoles
    variable = mean_en
    block = 0
    offset = 15.0
  []
  # [mean_en_advection_stabilization]
  #   type = EFieldArtDiff
  #   variable = mean_en
  #   potential = potential
  #   block = 0
  # []
[]

[Variables]
  [potential]
  []
  [em]
    block = 0
  []
  [emliq]
    block = 1
    # scaling = 1e-5
  []

  [Arp]
    block = 0
  []

  [mean_en]
    block = 0
    # scaling = 1e-1
  []

  [OHm]
    block = 1
    # scaling = 1e-5
  []
[]

[AuxVariables]
  [e_temp]
    block = 0
    order = CONSTANT
    family = MONOMIAL
  []
  #[x]
  #  order = CONSTANT
  #  family = MONOMIAL
  #[]
  [position0]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [position1]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [x_node]
  []
  [rho]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [rholiq]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  #[em_lin]
  #  order = CONSTANT
  #  family = MONOMIAL
  #  block = 0
  #[]
  [em_density]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  #[emliq_lin]
  #  order = CONSTANT
  #  family = MONOMIAL
  #  block = 1
  #[]
  [emliq_density]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  #[Arp_lin]
  #  order = CONSTANT
  #  family = MONOMIAL
  #  block = 0
  #[]
  [Arp_density]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  #[OHm_lin]
  #  block = 1
  #  order = CONSTANT
  #  family = MONOMIAL
  #[]
  [OHm_density]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  #[Efield]
  #  order = CONSTANT
  #  family = MONOMIAL
  #[]
  [EFieldx0]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [EFieldx1]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [Current_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [Current_emliq]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [Current_Arp]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [Current_OHm]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [tot_gas_current]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [tot_liq_current]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [tot_flux_OHm]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [EFieldAdvAux_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [DiffusiveFlux_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [EFieldAdvAux_emliq]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [DiffusiveFlux_emliq]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [PowerDep_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [PowerDep_Arp]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [ProcRate_el]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [ProcRate_ex]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [ProcRate_iz]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
[]

[AuxKernels]
  [PowerDep_em]
    type = ADPowerDep
    density_log = em
    potential = potential
    art_diff = false
    potential_units = kV
    variable = PowerDep_em
    position_units = ${dom0Scale}
    block = 0
  []
  [PowerDep_Arp]
    type = ADPowerDep
    density_log = Arp
    potential = potential
    art_diff = false
    potential_units = kV
    variable = PowerDep_Arp
    position_units = ${dom0Scale}
    block = 0
  []
  [ProcRate_el]
    type = ADProcRate
    em = em
    potential = potential
    proc = el
    variable = ProcRate_el
    position_units = ${dom0Scale}
    block = 0
  []
  [ProcRate_ex]
    type = ADProcRate
    em = em
    potential = potential
    proc = ex
    variable = ProcRate_ex
    position_units = ${dom0Scale}
    block = 0
  []
  [ProcRate_iz]
    type = ADProcRate
    em = em
    potential = potential
    proc = iz
    variable = ProcRate_iz
    position_units = ${dom0Scale}
    block = 0
  []
  [e_temp]
    type = ElectronTemperature
    variable = e_temp
    electron_density = em
    mean_en = mean_en
    block = 0
  []
  [x_g]
    type = Position
    variable = position0
    position_units = ${dom0Scale}
    block = 0
  []
  [x_l]
    type = Position
    variable = position1
    position_units = ${dom1Scale}
    block = 1
  []
  [x_ng]
    type = Position
    variable = x_node
    position_units = ${dom0Scale}
    block = 0
  []
  [x_nl]
    type = Position
    variable = x_node
    position_units = ${dom1Scale}
    block = 1
  []
  [rho]
    type = ParsedAux
    variable = rho
    coupled_variables = 'em_density Arp_density'
    expression = 'Arp_density - em_density'
    execute_on = 'timestep_end'
    block = 0
  []
  [rholiq]
    type = ParsedAux
    variable = rholiq
    coupled_variables = 'emliq_density OHm_density'
    expression = '-emliq_density - OHm_density'
    execute_on = 'timestep_end'
    block = 1
  []
  [tot_gas_current]
    type = ParsedAux
    variable = tot_gas_current
    coupled_variables = 'Current_em Current_Arp'
    expression = 'Current_em + Current_Arp'
    execute_on = 'timestep_end'
    block = 0
  []
  [tot_liq_current]
    type = ParsedAux
    variable = tot_liq_current
    coupled_variables = 'Current_emliq Current_OHm' # Current_H3Op Current_OHm'
    expression = 'Current_emliq + Current_OHm' # + Current_H3Op + Current_OHm'
    execute_on = 'timestep_end'
    block = 1
  []
  [em_lin]
    type = DensityMoles
    variable = em_density
    density_log = em
    block = 0
  []
  [emliq_lin]
    type = DensityMoles
    variable = emliq_density
    density_log = emliq
    block = 1
  []
  [Arp_lin]
    type = DensityMoles
    variable = Arp_density
    density_log = Arp
    block = 0
  []
  [OHm_lin]
    type = DensityMoles
    variable = OHm_density
    density_log = OHm
    block = 1
  []
  [Efield_g]
    type = Efield
    component = 0
    potential = potential
    variable = EFieldx0
    position_units = ${dom0Scale}
    block = 0
  []
  [Efield_l]
    type = Efield
    component = 0
    potential = potential
    variable = EFieldx1
    position_units = ${dom1Scale}
    block = 1
  []
  [Current_em]
    type = ADCurrent
    potential = potential
    density_log = em
    variable = Current_em
    art_diff = false
    block = 0
    position_units = ${dom0Scale}
  []
  [Current_emliq]
    type = ADCurrent
    potential = potential
    density_log = emliq
    variable = Current_emliq
    art_diff = false
    block = 1
    position_units = ${dom1Scale}
  []
  [Current_Arp]
    type = ADCurrent
    potential = potential
    density_log = Arp
    variable = Current_Arp
    art_diff = false
    block = 0
    position_units = ${dom0Scale}
  []
  [Current_OHm]
    block = 1
    type = ADCurrent
    potential = potential
    density_log = OHm
    variable = Current_OHm
    art_diff = false
    position_units = ${dom1Scale}
  []
  [tot_flux_OHm]
    block = 1
    type = ADTotalFlux
    potential = potential
    density_log = OHm
    variable = tot_flux_OHm
  []
  [EFieldAdvAux_em]
    type = ADEFieldAdvAux
    potential = potential
    density_log = em
    variable = EFieldAdvAux_em
    block = 0
    position_units = ${dom0Scale}
  []
  [DiffusiveFlux_em]
    type = ADDiffusiveFlux
    density_log = em
    variable = DiffusiveFlux_em
    block = 0
    position_units = ${dom0Scale}
  []
  [EFieldAdvAux_emliq]
    type = ADEFieldAdvAux
    potential = potential
    density_log = emliq
    variable = EFieldAdvAux_emliq
    block = 1
    position_units = ${dom1Scale}
  []
  [DiffusiveFlux_emliq]
    type = ADDiffusiveFlux
    density_log = emliq
    variable = DiffusiveFlux_emliq
    block = 1
    position_units = ${dom1Scale}
  []
[]

[InterfaceKernels]
  [em_advection]
    type = InterfaceAdvection
    mean_en_neighbor = mean_en
    potential_neighbor = potential
    neighbor_var = em
    variable = emliq
    boundary = master1_interface
    position_units = ${dom1Scale}
    neighbor_position_units = ${dom0Scale}
  []
  [em_diffusion]
    type = InterfaceLogDiffusionElectrons
    mean_en_neighbor = mean_en
    neighbor_var = em
    variable = emliq
    boundary = master1_interface
    position_units = ${dom1Scale}
    neighbor_position_units = ${dom0Scale}
  []
[]

[BCs]
  [mean_en_physical_right]
    type = HagelaarEnergyBC
    variable = mean_en
    boundary = 'master0_interface'
    potential = potential
    electrons = em
    r = 0.99
    position_units = ${dom0Scale}
  []
  [mean_en_physical_left]
    type = HagelaarEnergyBC
    variable = mean_en
    boundary = 'left'
    potential = potential
    electrons = em
    r = 0
    position_units = ${dom0Scale}
  []
  [secondary_energy_left]
    type = SecondaryElectronEnergyBC
    variable = mean_en
    boundary = 'left'
    potential = potential
    electrons = em
    ions = 'Arp'
    r = 0
    emission_coeffs = 0.05
    secondary_electron_energy = 3
    position_units = ${dom0Scale}
  []

  [potential_left]
    type = NeumannCircuitVoltageMoles_KV
    variable = potential
    boundary = left
    function = potential_bc_func
    ions = Arp
    data_provider = data_provider
    electrons = em
    electron_energy = mean_en
    r = 0
    position_units = ${dom0Scale}
    emission_coeffs = 0.05
  []
  [potential_dirichlet_right]
    type = DirichletBC
    variable = potential
    boundary = right
    value = 0
  []
  [em_physical_right]
    type = HagelaarElectronBC
    variable = em
    boundary = 'master0_interface'
    potential = potential
    electron_energy = mean_en
    r = 0.99
    position_units = ${dom0Scale}
  []
  [Arp_physical_right_diffusion]
    type = HagelaarIonDiffusionBC
    variable = Arp
    boundary = 'master0_interface'
    r = 0
    position_units = ${dom0Scale}
  []
  [Arp_physical_right_advection]
    type = HagelaarIonAdvectionBC
    variable = Arp
    boundary = 'master0_interface'
    potential = potential
    r = 0
    position_units = ${dom0Scale}
  []

  [em_physical_left]
    type = HagelaarElectronBC
    variable = em
    boundary = 'left'
    potential = potential
    electron_energy = mean_en
    r = 0
    position_units = ${dom0Scale}
  []
  [sec_electrons_left]
    type = SecondaryElectronBC
    variable = em
    boundary = 'left'
    potential = potential
    ions = Arp
    electron_energy = mean_en
    r = 0
    position_units = ${dom0Scale}
    emission_coeffs = 0.05
  []
  [Arp_physical_left_diffusion]
    type = HagelaarIonDiffusionBC
    variable = Arp
    boundary = 'left'
    r = 0
    position_units = ${dom0Scale}
  []
  [Arp_physical_left_advection]
    type = HagelaarIonAdvectionBC
    variable = Arp
    boundary = 'left'
    potential = potential
    r = 0
    position_units = ${dom0Scale}
  []

  [emliq_right]
    type = DCIonBC
    variable = emliq
    boundary = right
    potential = potential
    position_units = ${dom1Scale}
  []
  [OHm_physical]
    type = DCIonBC
    variable = OHm
    boundary = 'right'
    potential = potential
    position_units = ${dom1Scale}
  []
[]

[ICs]
  [em_ic]
    type = ConstantIC
    variable = em
    value = -21
    block = 0
  []
  [emliq_ic]
    type = ConstantIC
    variable = emliq
    value = -21
    block = 1
  []
  [Arp_ic]
    type = ConstantIC
    variable = Arp
    value = -21
    block = 0
  []
  [mean_en_ic]
    type = ConstantIC
    variable = mean_en
    value = -20
    block = 0
  []
  [OHm_ic]
    type = ConstantIC
    variable = OHm
    value = -15.6
    block = 1
  []
  [potential_ic]
    type = FunctionIC
    variable = potential
    function = potential_ic_func
  []
[]

[Functions]
  [potential_bc_func]
    type = ParsedFunction
    # expression = '1.25*tanh(1e6*t)'
    expression = -1.25
  []
  [potential_ic_func]
    type = ParsedFunction
    expression = '-1.25 * (1.0001e-3 - x)'
  []
[]

[Materials]
  [gas_block]
    type = Gas
    interp_trans_coeffs = true
    interp_elastic_coeff = true
    ramp_trans_coeffs = false
    em = em
    potential = potential
    ip = Arp
    mean_en = mean_en
    user_se_coeff = 0.05
    block = 0
    property_tables_file = td_argon_mean_en.txt
  []
  [water_block]
    type = Water
    block = 1
    potential = potential
  []
[]

(tutorial/tutorial05-PlasmaWaterInterface/DC_argon-With-Water.i)

#This tutorial is of an argon DC discharge over water at
#different reflection coefficients. Reflection coefficients
#for electron and mean energy at the water are on lines 297 and 315.

#Scaling for block 0 (Plasma)
dom0Scale = 1e-3
#caling for block 1 (Water)
dom1Scale = 1e-7

[GlobalParams]
  # off set for log stabilization, prevents log(0)
  offset = 20
  #Scales the potential by V or kV
  potential_units = kV
  #Converts density from #/m^3 to moles/m^3
  use_moles = true
[]

[Mesh]
  #Mesh is defined by a Gmsh file
  [file]
    type = FileMeshGenerator
    file = 'plasmaliquid.msh'
  []

  [interface]
    # interface allows for one directional fluxes
    # interface going from air to water
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = '0' # block where the flux is coming from
    paired_block = '1' # block where the flux is going to
    new_boundary = 'master0_interface'
    input = file # first input is where the msh is coming and it is then named air_to_water
  []
  [interface_again]
    # interface going from water to air
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = '1' # block where the flux is coming from
    paired_block = '0' # block where the flux is going to
    new_boundary = 'master1_interface'
    input = interface
  []
  #Renames all sides with the specified normal
  #For 1D, this is used to rename the end points of the mesh
  [left]
    type = SideSetsFromNormalsGenerator
    normals = '-1 0 0'
    new_boundary = 'left'
    input = interface_again
  []
  [right]
    type = SideSetsFromNormalsGenerator
    normals = '1 0 0'
    new_boundary = 'right'
    input = left
  []
[]

#Defines the problem type, such as FE, eigen value problem, etc.
[Problem]
  type = FEProblem
[]

#The potential needs to be defined outside of the Action,
#since it is present in both Blocks
#Declare any variable that is present in mulitple blocks in the variables section
[Variables]
  [potential]
  []
[]

[DriftDiffusionAction]
  [Plasma]
    #User define name for electrons (usually 'em')
    electrons = em
    #User define name for ions
    charged_particle = Arp
    #User define name for potential (usually 'potential')
    potential = potential
    #Set False becuase both areas use the same potential
    Is_potential_unique = false
    #User define name for the electron mean energy density (usually 'mean_en')
    mean_energy = mean_en
    #Helps prevent the log(0)
    using_offset = true #helps prevent the log(0)
    #The position scaling for the mesh, define at top of input file
    position_units = ${dom0Scale}
    #Name of material block for plasma
    block = 0
    #Additional outputs, such as ElectronTemperature, Current, and EField.
    Additional_Outputs = 'ElectronTemperature Current EField'
  []
  # treats water as a dense plasma
  [Water]
    charged_particle = 'emliq OHm'
    potential = potential
    Is_potential_unique = false
    using_offset = true
    position_units = ${dom1Scale}
    #Name of material block for water
    block = 1
    Additional_Outputs = 'Current EField'
  []
[]

[Reactions]
  [Argon]
    #Name of reactant species that are variables
    species = 'em Arp'
    #Name of reactant species that are auxvariables
    aux_species = 'Ar'
    #Type of coefficient (rate or townsend)
    reaction_coefficient_format = 'townsend'
    #Name of background gas
    gas_species = 'Ar'
    #Name of the electron mean energy density (usually 'mean_en')
    electron_energy = 'mean_en'
    #Name of the electrons (usually 'em')
    electron_density = 'em'
    #Defines if electrons are tracked
    include_electrons = true
    #Defines directory holding rate text files
    file_location = 'townsend_coefficients'
    #Name of name for potential (usually 'potential')
    potential = 'potential'
    #Defines if log form is used (true for Zapdos)
    use_log = true
    #Defines if automatic differentiation is used (true for Zapdos)
    use_ad = true
    #The position scaling for the mesh, define at top of input file
    position_units = ${dom0Scale}
    #Name of material block for plasma
    block = 0
    #Inputs of the plasma chemsity
    #e.g. Reaction : Constant or EEDF dependent [Threshold Energy] (Text file name)
    #     em + Ar -> em + Ar*        : EEDF [-11.56] (reaction1)
    reactions = 'em + Ar -> em + Ar               : EEDF [elastic] (reaction1)
                 em + Ar -> em + Ar*              : EEDF [-11.5] (reaction2)
                 em + Ar -> em + em + Arp         : EEDF [-15.76] (reaction3)'
  []

  [Water]
    species = 'emliq OHm'
    reaction_coefficient_format = 'rate'
    use_log = true
    use_ad = true
    aux_species = 'H2O'
    block = 1
    reactions = 'emliq -> H + OHm : 1064
                 emliq + emliq -> H2 + OHm + OHm : 3.136e8'
  []
[]

[AuxVariables]
  #Background fluid in water
  [H2O]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 10.92252
    block = 1
  []
  #Background gas in plasma
  [Ar]
    block = 0
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 3.70109
  []
[]

[InterfaceKernels]
  # interface conditions allow one this to go another
  # These account for electrons going into the water
  # if you want to account for electrons going out of the water into the air then you would need there to be interface conditions for master interface 0
  [em_advection]
    type = InterfaceAdvection
    mean_en_neighbor = mean_en
    potential_neighbor = potential
    neighbor_var = em

    #the main variable being affected. The em is going into the water so em -> emliq
    # this affects emliq
    variable = emliq
    boundary = master1_interface
    position_units = ${dom1Scale}
    neighbor_position_units = ${dom0Scale}
  []
  [em_diffusion]
    type = InterfaceLogDiffusionElectrons
    mean_en_neighbor = mean_en
    neighbor_var = em
    variable = emliq
    boundary = master1_interface
    position_units = ${dom1Scale}
    neighbor_position_units = ${dom0Scale}
  []
[]

#Electrode data needed for 1D BC of 'NeumannCircuitVoltageMoles_KV'
[UserObjects]
  [data_provider]
    type = ProvideMobility
    electrode_area = 5.02e-7 # Formerly 3.14e-6
    ballast_resist = 1e6
    e = 1.6e-19
  []
[]

[BCs]
  #DC BC on the air electrode
  [potential_left]
    type = NeumannCircuitVoltageMoles_KV
    variable = potential
    boundary = left
    function = potential_bc_func
    ions = Arp
    data_provider = data_provider
    electrons = em
    electron_energy = mean_en
    r = 0
    position_units = ${dom0Scale}
    emission_coeffs = 0.05
  []
  #Ground electrode under the water
  [potential_dirichlet_right]
    type = DirichletBC
    variable = potential
    boundary = right
    value = 0
  []

  #Electrons on the powered electrode
  [em_physical_left]
    type = HagelaarElectronBC
    variable = em
    boundary = 'left'
    potential = potential
    electron_energy = mean_en
    r = 0
    position_units = ${dom0Scale}
  []
  [sec_electrons_left]
    type = SecondaryElectronBC
    variable = em
    boundary = 'left'
    potential = potential
    ions = Arp
    electron_energy = mean_en
    r = 0
    position_units = ${dom0Scale}
    emission_coeffs = 0.05
  []

  #Mean electron energy on the powered electrode
  [mean_en_physical_left]
    type = HagelaarEnergyBC
    variable = mean_en
    boundary = 'left'
    potential = potential
    electrons = em
    r = 0
    position_units = ${dom0Scale}
  []
  [secondary_energy_left]
    type = SecondaryElectronEnergyBC
    variable = mean_en
    boundary = 'left'
    potential = potential
    electrons = em
    ions = 'Arp'
    r = 0
    position_units = ${dom0Scale}
    emission_coeffs = 0.05
    secondary_electron_energy = 3
  []

  #Argon ions on the powered electrode
  [Arp_physical_left_diffusion]
    type = HagelaarIonDiffusionBC
    variable = Arp
    boundary = 'left'
    r = 0
    position_units = ${dom0Scale}
  []
  [Arp_physical_left_advection]
    type = HagelaarIonAdvectionBC
    variable = Arp
    boundary = 'left'
    potential = potential
    r = 0
    position_units = ${dom0Scale}
  []

  #Electrons on the plasma-liquid interface
  [em_physical_right]
    type = HagelaarElectronBC
    variable = em
    boundary = 'master0_interface'
    potential = potential
    electron_energy = mean_en
    r = 0.00
    position_units = ${dom0Scale}
  []
  [Arp_physical_right_diffusion]
    type = HagelaarIonDiffusionBC
    variable = Arp
    boundary = 'master0_interface'
    r = 0
    position_units = ${dom0Scale}
  []
  [mean_en_physical_right]
    #Mean electron energy on the plasma-liquid interface
    type = HagelaarEnergyBC
    variable = mean_en
    boundary = 'master0_interface'
    potential = potential
    electrons = em
    r = 0.00
    position_units = ${dom0Scale}
  []

  #Argon ions on the plasma-liquid interface
  [Arp_physical_right_advection]
    type = HagelaarIonAdvectionBC
    variable = Arp
    boundary = 'master0_interface'
    potential = potential
    r = 0
    position_units = ${dom0Scale}
  []

  #Electrons on the electrode
  [emliq_right]
    type = DCIonBC
    variable = emliq
    boundary = 'right'
    potential = potential
    position_units = ${dom1Scale}
  []
  #OH- on the ground electrode
  [OHm_physical]
    type = DCIonBC
    variable = OHm
    boundary = 'right'
    potential = potential
    position_units = ${dom1Scale}
  []
[]

#Initial conditions for variables.
#If left undefine, the IC is zero
[ICs]
  [em_ic]
    type = ConstantIC
    variable = em
    value = -21
    block = 0
  []
  [emliq_ic]
    type = ConstantIC
    variable = emliq
    value = -21
    block = 1
  []
  [Arp_ic]
    type = ConstantIC
    variable = Arp
    value = -21
    block = 0
  []
  [mean_en_ic]
    type = ConstantIC
    variable = mean_en
    value = -20
    block = 0
  []
  [OHm_ic]
    type = ConstantIC
    variable = OHm
    value = -15.6
    block = 1
  []
  [potential_ic]
    type = FunctionIC
    variable = potential
    function = potential_ic_func
  []
[]

#Define function used throughout the input file (e.g. BCs and ICs)
[Functions]
  [potential_bc_func]
    type = ParsedFunction
    # expression = '1.25*tanh(1e6*t)'
    expression = -1.25
  []
  [potential_ic_func]
    type = ParsedFunction
    expression = '-1.25 * (1.0001e-3 - x)'
  []
[]

[Materials]
  #The material properties for electrons and ions in water
  [water_block]
    type = Water
    block = 1
    potential = potential
  []

  #The material properties for electrons in plasma
  #Also hold universal constant, such as Avogadro's number, elementary charge, etc.
  [electrons_in_plasma]
    type = GasElectronMoments
    interp_trans_coeffs = true
    interp_elastic_coeff = true
    ramp_trans_coeffs = false
    user_p_gas = 101325
    user_se_coeff = 0.05
    em = em
    potential = potential
    mean_en = mean_en
    block = 0
    property_tables_file = 'townsend_coefficients/moments.txt'
  []

  #Sets the pressure and temperature in the water
  [water_block1]
    type = GenericConstantMaterial
    block = 1
    prop_names = 'T_gas p_gas'
    prop_values = '300 1.01e5'
  []

  #The material properties of the argon ion
  [gas_species_0]
    type = ADHeavySpecies
    heavy_species_name = Arp
    heavy_species_mass = 6.64e-26
    heavy_species_charge = 1.0
    block = 0
  []

  #The material properties of the background argon gas
  [gas_species_2]
    type = ADHeavySpecies
    heavy_species_name = Ar
    heavy_species_mass = 6.64e-26
    heavy_species_charge = 0.0
    block = 0
  []
[]

#Preconditioning options
#Learn more at: https://mooseframework.inl.gov/syntax/Preconditioning/index.html
[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

#How to execute the problem.
#Defines type of solve (such as steady or transient),
# solve type (Newton, PJFNK, etc.) and tolerances
[Executioner]
  type = Transient
  end_time = 1e-1
  petsc_options = '-snes_converged_reason -snes_linesearch_monitor'
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount'
  petsc_options_value = 'lu NONZERO 1.e-10'
  nl_rel_tol = 1e-4
  nl_abs_tol = 7.6e-5
  dtmin = 1e-15
  l_max_its = 20
  [TimeSteppers]
    [Adaptive]
      type = IterationAdaptiveDT
      cutback_factor = 0.4
      dt = 1e-11
      growth_factor = 1.2
      optimal_iterations = 30
    []
  []
[]

#Defines the output type of the file (multiple output files can be define per run)
[Outputs]
  perf_graph = true
  [out]
    type = Exodus
    execute_on = 'final'
  []
[]

(test/tests/1d_dc/mean_en.i)

dom0Scale = 1e-3
dom1Scale = 1e-7

[GlobalParams]
  offset = 20
  # offset = 0
  potential_units = kV
  use_moles = true
  # potential_units = V
[]

[Mesh]
  [file]
    type = FileMeshGenerator
    file = 'liquidNew.msh'
  []
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = '0'
    paired_block = '1'
    new_boundary = 'master0_interface'
    input = file
  []
  [interface_again]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = '1'
    paired_block = '0'
    new_boundary = 'master1_interface'
    input = interface
  []
  [left]
    type = SideSetsFromNormalsGenerator
    normals = '-1 0 0'
    new_boundary = 'left'
    input = interface_again
  []
  [right]
    type = SideSetsFromNormalsGenerator
    normals = '1 0 0'
    new_boundary = 'right'
    input = left
  []
[]

[Problem]
  type = FEProblem
  # kernel_coverage_check = false
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
    ksp_norm = none
  []
[]

[Executioner]
  type = Transient
  end_time = 1e-1
  petsc_options = '-snes_converged_reason -snes_linesearch_monitor'
  #petsc_options = '-snes_converged_reason -snes_linesearch_monitor -snes_test_jacobian'
  # petsc_options = '-snes_test_display'
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount -ksp_type -snes_linesearch_minlambda'
  petsc_options_value = 'lu NONZERO 1.e-10 preonly 1e-3'
  # petsc_options_iname = '-pc_type -sub_pc_type'
  # petsc_options_value = 'asm lu'
  # petsc_options_iname = '-snes_type'
  # petsc_options_value = 'test'
  nl_rel_tol = 1e-4
  nl_abs_tol = 7.6e-5
  dtmin = 1e-12
  l_max_its = 20
  [TimeSteppers]
    [Adaptive]
      type = IterationAdaptiveDT
      cutback_factor = 0.4
      dt = 1e-11
      # dt = 1.1
      growth_factor = 1.2
      optimal_iterations = 15
    []
  []
[]

[Outputs]
  perf_graph = true
  # print_linear_residuals = false
  [out]
    type = Exodus
    execute_on = 'final'
  []
  [dof_map]
    type = DOFMap
  []
[]

[Debug]
  show_var_residual_norms = true
[]

[UserObjects]
  [data_provider]
    type = ProvideMobility
    electrode_area = 5.02e-7 # Formerly 3.14e-6
    ballast_resist = 1e6
    e = 1.6e-19
  []
[]

[Kernels]
  [em_time_deriv]
    type = ElectronTimeDerivative
    variable = em
    block = 0
  []
  [em_advection]
    type = EFieldAdvection
    variable = em
    potential = potential
    block = 0
    position_units = ${dom0Scale}
  []
  [em_diffusion]
    type = CoeffDiffusion
    variable = em
    block = 0
    position_units = ${dom0Scale}
  []
  [em_ionization]
    type = ElectronsFromIonization
    variable = em
    potential = potential
    mean_en = mean_en
    em = em
    block = 0
    position_units = ${dom0Scale}
  []
  [em_log_stabilization]
    type = LogStabilizationMoles
    variable = em
    block = 0
  []

  [emliq_time_deriv]
    type = ElectronTimeDerivative
    variable = emliq
    block = 1
  []
  [emliq_advection]
    type = EFieldAdvection
    variable = emliq
    potential = potential
    block = 1
    position_units = ${dom1Scale}
  []
  [emliq_diffusion]
    type = CoeffDiffusion
    variable = emliq
    block = 1
    position_units = ${dom1Scale}
  []
  [emliq_reactant_first_order_rxn]
    type = ReactantFirstOrderRxn
    variable = emliq
    block = 1
  []
  [emliq_water_bi_sink]
    type = ReactantAARxn
    variable = emliq
    block = 1
  []
  [emliq_log_stabilization]
    type = LogStabilizationMoles
    variable = emliq
    block = 1
  []

  [potential_diffusion_dom1]
    type = CoeffDiffusionLin
    variable = potential
    block = 0
    position_units = ${dom0Scale}
  []
  [potential_diffusion_dom2]
    type = CoeffDiffusionLin
    variable = potential
    block = 1
    position_units = ${dom1Scale}
  []
  [Arp_charge_source]
    type = ChargeSourceMoles_KV
    variable = potential
    charged = Arp
    block = 0
  []
  [em_charge_source]
    type = ChargeSourceMoles_KV
    variable = potential
    charged = em
    block = 0
  []
  [emliq_charge_source]
    type = ChargeSourceMoles_KV
    variable = potential
    charged = emliq
    block = 1
  []
  [OHm_charge_source]
    type = ChargeSourceMoles_KV
    variable = potential
    charged = OHm
    block = 1
  []

  [Arp_time_deriv]
    type = ElectronTimeDerivative
    variable = Arp
    block = 0
  []
  [Arp_advection]
    type = EFieldAdvection
    variable = Arp
    potential = potential
    position_units = ${dom0Scale}
    block = 0
  []
  [Arp_diffusion]
    type = CoeffDiffusion
    variable = Arp
    block = 0
    position_units = ${dom0Scale}
  []
  [Arp_ionization]
    type = IonsFromIonization
    variable = Arp
    potential = potential
    em = em
    mean_en = mean_en
    block = 0
    position_units = ${dom0Scale}
  []
  [Arp_log_stabilization]
    type = LogStabilizationMoles
    variable = Arp
    block = 0
  []

  [OHm_time_deriv]
    type = ElectronTimeDerivative
    variable = OHm
    block = 1
  []
  [OHm_advection]
    type = EFieldAdvection
    variable = OHm
    potential = potential
    block = 1
    position_units = ${dom1Scale}
  []
  [OHm_diffusion]
    type = CoeffDiffusion
    variable = OHm
    block = 1
    position_units = ${dom1Scale}
  []
  [OHm_log_stabilization]
    type = LogStabilizationMoles
    variable = OHm
    block = 1
  []
  [OHm_product_first_order_rxn]
    type = ProductFirstOrderRxn
    variable = OHm
    v = emliq
    block = 1
  []
  [OHm_product_aabb_rxn]
    type = ProductAABBRxn
    variable = OHm
    v = emliq
    block = 1
  []

  [mean_en_time_deriv]
    type = ElectronTimeDerivative
    variable = mean_en
    block = 0
  []
  [mean_en_advection]
    type = EFieldAdvection
    variable = mean_en
    potential = potential
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_diffusion]
    type = CoeffDiffusion
    variable = mean_en
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_joule_heating]
    type = JouleHeating
    variable = mean_en
    potential = potential
    em = em
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_ionization]
    type = ElectronEnergyLossFromIonization
    variable = mean_en
    potential = potential
    em = em
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_elastic]
    type = ElectronEnergyLossFromElastic
    variable = mean_en
    potential = potential
    em = em
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_excitation]
    type = ElectronEnergyLossFromExcitation
    variable = mean_en
    potential = potential
    em = em
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_log_stabilization]
    type = LogStabilizationMoles
    variable = mean_en
    block = 0
    offset = 15
  []
[]

[Variables]
  [potential]
  []
  [em]
    block = 0
  []
  [emliq]
    block = 1
    # scaling = 1e-5
  []

  [Arp]
    block = 0
  []

  [mean_en]
    block = 0
    # scaling = 1e-1
  []

  [OHm]
    block = 1
    # scaling = 1e-5
  []
[]

[AuxVariables]
  [e_temp]
    block = 0
    order = CONSTANT
    family = MONOMIAL
  []
  [x]
    order = CONSTANT
    family = MONOMIAL
  []
  [x_node]
  []
  [rho]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [rholiq]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [em_lin]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [emliq_lin]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [Arp_lin]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [OHm_lin]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [Efield]
    order = CONSTANT
    family = MONOMIAL
  []
  [Current_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [Current_emliq]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [Current_Arp]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [Current_OHm]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [tot_gas_current]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [tot_liq_current]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [tot_flux_OHm]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [EFieldAdvAux_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [DiffusiveFlux_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [EFieldAdvAux_emliq]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [DiffusiveFlux_emliq]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [PowerDep_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [PowerDep_Arp]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [ProcRate_el]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [ProcRate_ex]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [ProcRate_iz]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
[]

[AuxKernels]
  [PowerDep_em]
    type = ADPowerDep
    density_log = em
    potential = potential
    art_diff = false
    potential_units = kV
    variable = PowerDep_em
    position_units = ${dom0Scale}
    block = 0
  []
  [PowerDep_Arp]
    type = ADPowerDep
    density_log = Arp
    potential = potential
    art_diff = false
    potential_units = kV
    variable = PowerDep_Arp
    position_units = ${dom0Scale}
    block = 0
  []
  [ProcRate_el]
    type = ADProcRate
    em = em
    potential = potential
    proc = el
    variable = ProcRate_el
    position_units = ${dom0Scale}
    block = 0
  []
  [ProcRate_ex]
    type = ADProcRate
    em = em
    potential = potential
    proc = ex
    variable = ProcRate_ex
    position_units = ${dom0Scale}
    block = 0
  []
  [ProcRate_iz]
    type = ADProcRate
    em = em
    potential = potential
    proc = iz
    variable = ProcRate_iz
    position_units = ${dom0Scale}
    block = 0
  []
  [e_temp]
    type = ElectronTemperature
    variable = e_temp
    electron_density = em
    mean_en = mean_en
    block = 0
  []
  [x_g]
    type = Position
    variable = x
    position_units = ${dom0Scale}
    block = 0
  []
  [x_l]
    type = Position
    variable = x
    position_units = ${dom1Scale}
    block = 1
  []
  [x_ng]
    type = Position
    variable = x_node
    position_units = ${dom0Scale}
    block = 0
  []
  [x_nl]
    type = Position
    variable = x_node
    position_units = ${dom1Scale}
    block = 1
  []
  [rho]
    type = ParsedAux
    variable = rho
    coupled_variables = 'em_lin Arp_lin'
    expression = 'Arp_lin - em_lin'
    execute_on = 'timestep_end'
    block = 0
  []
  [rholiq]
    type = ParsedAux
    variable = rholiq
    coupled_variables = 'emliq_lin OHm_lin' # H3Op_lin OHm_lin'
    expression = '-emliq_lin - OHm_lin' # 'H3Op_lin - em_lin - OHm_lin'
    execute_on = 'timestep_end'
    block = 1
  []
  [tot_gas_current]
    type = ParsedAux
    variable = tot_gas_current
    coupled_variables = 'Current_em Current_Arp'
    expression = 'Current_em + Current_Arp'
    execute_on = 'timestep_end'
    block = 0
  []
  [tot_liq_current]
    type = ParsedAux
    variable = tot_liq_current
    coupled_variables = 'Current_emliq Current_OHm' # Current_H3Op Current_OHm'
    expression = 'Current_emliq + Current_OHm' # + Current_H3Op + Current_OHm'
    execute_on = 'timestep_end'
    block = 1
  []
  [em_lin]
    type = DensityMoles
    variable = em_lin
    density_log = em
    block = 0
  []
  [emliq_lin]
    type = DensityMoles
    variable = emliq_lin
    density_log = emliq
    block = 1
  []
  [Arp_lin]
    type = DensityMoles
    variable = Arp_lin
    density_log = Arp
    block = 0
  []
  [OHm_lin]
    type = DensityMoles
    variable = OHm_lin
    density_log = OHm
    block = 1
  []
  [Efield_g]
    type = Efield
    component = 0
    potential = potential
    variable = Efield
    position_units = ${dom0Scale}
    block = 0
  []
  [Efield_l]
    type = Efield
    component = 0
    potential = potential
    variable = Efield
    position_units = ${dom1Scale}
    block = 1
  []
  [Current_em]
    type = ADCurrent
    potential = potential
    density_log = em
    variable = Current_em
    art_diff = false
    block = 0
    position_units = ${dom0Scale}
  []
  [Current_emliq]
    type = ADCurrent
    potential = potential
    density_log = emliq
    variable = Current_emliq
    art_diff = false
    block = 1
    position_units = ${dom1Scale}
  []
  [Current_Arp]
    type = ADCurrent
    potential = potential
    density_log = Arp
    variable = Current_Arp
    art_diff = false
    block = 0
    position_units = ${dom0Scale}
  []
  [Current_OHm]
    block = 1
    type = ADCurrent
    potential = potential
    density_log = OHm
    variable = Current_OHm
    art_diff = false
    position_units = ${dom1Scale}
  []
  [tot_flux_OHm]
    block = 1
    type = ADTotalFlux
    potential = potential
    density_log = OHm
    variable = tot_flux_OHm
  []
  [EFieldAdvAux_em]
    type = ADEFieldAdvAux
    potential = potential
    density_log = em
    variable = EFieldAdvAux_em
    block = 0
    position_units = ${dom0Scale}
  []
  [DiffusiveFlux_em]
    type = ADDiffusiveFlux
    density_log = em
    variable = DiffusiveFlux_em
    block = 0
    position_units = ${dom0Scale}
  []
  [EFieldAdvAux_emliq]
    type = ADEFieldAdvAux
    potential = potential
    density_log = emliq
    variable = EFieldAdvAux_emliq
    block = 1
    position_units = ${dom1Scale}
  []
  [DiffusiveFlux_emliq]
    type = ADDiffusiveFlux
    density_log = emliq
    variable = DiffusiveFlux_emliq
    block = 1
    position_units = ${dom1Scale}
  []
[]

[InterfaceKernels]
  [em_advection]
    type = InterfaceAdvection
    mean_en_neighbor = mean_en
    potential_neighbor = potential
    neighbor_var = em
    variable = emliq
    boundary = master1_interface
    position_units = ${dom1Scale}
    neighbor_position_units = ${dom0Scale}
  []
  [em_diffusion]
    type = InterfaceLogDiffusionElectrons
    mean_en_neighbor = mean_en
    neighbor_var = em
    variable = emliq
    boundary = master1_interface
    position_units = ${dom1Scale}
    neighbor_position_units = ${dom0Scale}
  []
[]

[BCs]
  [mean_en_physical_right]
    type = HagelaarEnergyBC
    variable = mean_en
    boundary = 'master0_interface'
    potential = potential
    electrons = em
    r = 0.99
    position_units = ${dom0Scale}
  []
  [mean_en_physical_left]
    type = HagelaarEnergyBC
    variable = mean_en
    boundary = 'left'
    potential = potential
    electrons = em
    r = 0
    position_units = ${dom0Scale}
  []
  [secondary_energy_left]
    type = SecondaryElectronEnergyBC
    variable = mean_en
    boundary = 'left'
    potential = potential
    electrons = em
    ions = 'Arp'
    r = 0
    emission_coeffs = 0.05
    secondary_electron_energy = 3
    position_units = ${dom0Scale}
  []

  [potential_left]
    type = NeumannCircuitVoltageMoles_KV
    variable = potential
    boundary = left
    function = potential_bc_func
    ions = Arp
    data_provider = data_provider
    electrons = em
    electron_energy = mean_en
    r = 0
    position_units = ${dom0Scale}
    emission_coeffs = 0.05
  []
  [potential_dirichlet_right]
    type = DirichletBC
    variable = potential
    boundary = right
    value = 0
    preset = false
  []
  [em_physical_right]
    type = HagelaarElectronBC
    variable = em
    boundary = 'master0_interface'
    potential = potential
    electron_energy = mean_en
    r = 0.99
    position_units = ${dom0Scale}
  []
  [Arp_physical_right_diffusion]
    type = HagelaarIonDiffusionBC
    variable = Arp
    boundary = 'master0_interface'
    r = 0
    position_units = ${dom0Scale}
  []
  [Arp_physical_right_advection]
    type = HagelaarIonAdvectionBC
    variable = Arp
    boundary = 'master0_interface'
    potential = potential
    r = 0
    position_units = ${dom0Scale}
  []

  [em_physical_left]
    type = HagelaarElectronBC
    variable = em
    boundary = 'left'
    potential = potential
    electron_energy = mean_en
    r = 0
    position_units = ${dom0Scale}
  []
  [sec_electrons_left]
    type = SecondaryElectronBC
    variable = em
    boundary = 'left'
    potential = potential
    ions = Arp
    electron_energy = mean_en
    r = 0
    position_units = ${dom0Scale}
    emission_coeffs = 0.05
  []
  [Arp_physical_left_diffusion]
    type = HagelaarIonDiffusionBC
    variable = Arp
    boundary = 'left'
    r = 0
    position_units = ${dom0Scale}
  []
  [Arp_physical_left_advection]
    type = HagelaarIonAdvectionBC
    variable = Arp
    boundary = 'left'
    potential = potential
    r = 0
    position_units = ${dom0Scale}
  []

  [emliq_right]
    type = DCIonBC
    variable = emliq
    boundary = right
    potential = potential
    position_units = ${dom1Scale}
  []
  [OHm_physical]
    type = DCIonBC
    variable = OHm
    boundary = 'right'
    potential = potential
    position_units = ${dom1Scale}
  []
[]

[ICs]
  [em_ic]
    type = ConstantIC
    variable = em
    value = -21
    block = 0
  []
  [emliq_ic]
    type = ConstantIC
    variable = emliq
    value = -21
    block = 1
  []
  [Arp_ic]
    type = ConstantIC
    variable = Arp
    value = -21
    block = 0
  []
  [mean_en_ic]
    type = ConstantIC
    variable = mean_en
    value = -20
    block = 0
  []
  [OHm_ic]
    type = ConstantIC
    variable = OHm
    value = -15.6
    block = 1
  []
  [potential_ic]
    type = FunctionIC
    variable = potential
    function = potential_ic_func
  []
[]

[Functions]
  [potential_bc_func]
    type = ParsedFunction
    # expression = '1.25*tanh(1e6*t)'
    expression = -1.25
  []
  [potential_ic_func]
    type = ParsedFunction
    expression = '-1.25 * (1.0001e-3 - x)'
  []
[]

[Materials]
  [gas_block]
    type = Gas
    interp_trans_coeffs = true
    interp_elastic_coeff = true
    ramp_trans_coeffs = false
    em = em
    potential = potential
    ip = Arp
    mean_en = mean_en
    user_se_coeff = 0.05
    block = 0
    property_tables_file = td_argon_mean_en.txt
  []
  [water_block]
    type = Water
    block = 1
    potential = potential
  []
[]

(test/tests/DriftDiffusionAction/mean_en_actions.i)

dom0Scale = 1e-3
dom1Scale = 1e-7

[GlobalParams]
  offset = 20.0
  # offset = 0
  potential_units = kV
  use_moles = true
  # potential_units = V
[]

[Mesh]
  [file]
    type = FileMeshGenerator
    file = 'liquidNew.msh'
  []
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = '0'
    paired_block = '1'
    new_boundary = 'master0_interface'
    input = file
  []
  [interface_again]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = '1'
    paired_block = '0'
    new_boundary = 'master1_interface'
    input = interface
  []
  [left]
    type = SideSetsFromNormalsGenerator
    normals = '-1 0 0'
    new_boundary = 'left'
    input = interface_again
  []
  [right]
    type = SideSetsFromNormalsGenerator
    normals = '1 0 0'
    new_boundary = 'right'
    input = left
  []
[]

[Problem]
  type = FEProblem
  # kernel_coverage_check = false
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
    ksp_norm = none
  []
[]

[Executioner]
  type = Transient
  end_time = 1e-1
  petsc_options = '-snes_converged_reason -snes_linesearch_monitor'
  # petsc_options = '-snes_test_display'
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount -ksp_type -snes_linesearch_minlambda'
  petsc_options_value = 'lu NONZERO 1.e-10 preonly 1e-3'
  # petsc_options_iname = '-pc_type -sub_pc_type'
  # petsc_options_value = 'asm lu'
  # petsc_options_iname = '-snes_type'
  # petsc_options_value = 'test'
  nl_rel_tol = 1e-4
  nl_abs_tol = 7.6e-5
  dtmin = 1e-12
  l_max_its = 20
  [TimeSteppers]
    [Adaptive]
      type = IterationAdaptiveDT
      cutback_factor = 0.4
      dt = 1e-11
      # dt = 1.1
      growth_factor = 1.2
      optimal_iterations = 15
    []
  []
[]

[Outputs]
  file_base = out
  perf_graph = true
  [out]
    type = Exodus
    execute_on = 'final'
  []
[]

[Debug]
  show_var_residual_norms = true
[]

[UserObjects]
  [data_provider]
    type = ProvideMobility
    electrode_area = 5.02e-7 # Formerly 3.14e-6
    ballast_resist = 1e6
    e = 1.6e-19
    # electrode_area = 1.1
    # ballast_resist = 1.1
    # e = 1.1
  []
[]

#The potential needs to be defined outside of the Action,
#since it is present in both Blocks
[Variables]
  [potential]
  []
[]

#Action the supplies the drift-diffusion equations for both Blocks,
#This action also adds JouleHeating and the ChargeSourceMoles_KV Kernels
[DriftDiffusionAction]
  [Plasma]
    electrons = em
    charged_particle = Arp
    potential = potential
    Is_potential_unique = false
    mean_energy = mean_en
    using_offset = true
    position_units = ${dom0Scale}
    block = 0
    Additional_Outputs = 'ElectronTemperature Current EField'
  []
  [Water]
    electrons = emliq
    charged_particle = OHm
    potential = potential
    Is_potential_unique = false
    using_offset = true
    position_units = ${dom1Scale}
    block = 1
    Additional_Outputs = 'Current EField'
  []
[]

#The Kernels supply the sources terms
[Kernels]
  [em_ionization]
    type = ElectronsFromIonization
    variable = em
    potential = potential
    mean_en = mean_en
    em = em
    block = 0
    position_units = ${dom0Scale}
  []

  [emliq_reactant_first_order_rxn]
    type = ReactantFirstOrderRxn
    variable = emliq
    block = 1
  []
  [emliq_water_bi_sink]
    type = ReactantAARxn
    variable = emliq
    block = 1
  []

  [Arp_ionization]
    type = IonsFromIonization
    variable = Arp
    potential = potential
    em = em
    mean_en = mean_en
    block = 0
    position_units = ${dom0Scale}
  []

  [OHm_product_first_order_rxn]
    type = ProductFirstOrderRxn
    variable = OHm
    v = emliq
    block = 1
  []
  [OHm_product_aabb_rxn]
    type = ProductAABBRxn
    variable = OHm
    v = emliq
    block = 1
  []

  [mean_en_ionization]
    type = ElectronEnergyLossFromIonization
    variable = mean_en
    potential = potential
    em = em
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_elastic]
    type = ElectronEnergyLossFromElastic
    variable = mean_en
    potential = potential
    em = em
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_excitation]
    type = ElectronEnergyLossFromExcitation
    variable = mean_en
    potential = potential
    em = em
    block = 0
    position_units = ${dom0Scale}
  []
[]

[AuxVariables]
  [x_node]
  []
  [rho]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [rholiq]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [tot_gas_current]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [tot_liq_current]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [tot_flux_OHm]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [EFieldAdvAux_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [DiffusiveFlux_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [EFieldAdvAux_emliq]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [DiffusiveFlux_emliq]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [PowerDep_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [PowerDep_Arp]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [ProcRate_el]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [ProcRate_ex]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [ProcRate_iz]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
[]

[AuxKernels]
  [PowerDep_em]
    type = ADPowerDep
    density_log = em
    potential = potential
    art_diff = false
    potential_units = kV
    variable = PowerDep_em
    position_units = ${dom0Scale}
    block = 0
  []
  [PowerDep_Arp]
    type = ADPowerDep
    density_log = Arp
    potential = potential
    art_diff = false
    potential_units = kV
    variable = PowerDep_Arp
    position_units = ${dom0Scale}
    block = 0
  []
  [ProcRate_el]
    type = ADProcRate
    em = em
    potential = potential
    proc = el
    variable = ProcRate_el
    position_units = ${dom0Scale}
    block = 0
  []
  [ProcRate_ex]
    type = ADProcRate
    em = em
    potential = potential
    proc = ex
    variable = ProcRate_ex
    position_units = ${dom0Scale}
    block = 0
  []
  [ProcRate_iz]
    type = ADProcRate
    em = em
    potential = potential
    proc = iz
    variable = ProcRate_iz
    position_units = ${dom0Scale}
    block = 0
  []
  [x_ng]
    type = Position
    variable = x_node
    position_units = ${dom0Scale}
    block = 0
  []
  [x_nl]
    type = Position
    variable = x_node
    position_units = ${dom1Scale}
    block = 1
  []
  [rho]
    type = ParsedAux
    variable = rho
    coupled_variables = 'em_density Arp_density'
    expression = 'Arp_density - em_density'
    execute_on = 'timestep_end'
    block = 0
  []
  [rholiq]
    type = ParsedAux
    variable = rholiq
    coupled_variables = 'emliq_density OHm_density'
    expression = '-emliq_density - OHm_density'
    execute_on = 'timestep_end'
    block = 1
  []
  [tot_gas_current]
    type = ParsedAux
    variable = tot_gas_current
    coupled_variables = 'Current_em Current_Arp'
    expression = 'Current_em + Current_Arp'
    execute_on = 'timestep_end'
    block = 0
  []
  [tot_liq_current]
    type = ParsedAux
    variable = tot_liq_current
    coupled_variables = 'Current_emliq Current_OHm' # Current_H3Op Current_OHm'
    expression = 'Current_emliq + Current_OHm' # + Current_H3Op + Current_OHm'
    execute_on = 'timestep_end'
    block = 1
  []
  [tot_flux_OHm]
    block = 1
    type = ADTotalFlux
    potential = potential
    density_log = OHm
    variable = tot_flux_OHm
  []
  [EFieldAdvAux_em]
    type = ADEFieldAdvAux
    potential = potential
    density_log = em
    variable = EFieldAdvAux_em
    block = 0
    position_units = ${dom0Scale}
  []
  [DiffusiveFlux_em]
    type = ADDiffusiveFlux
    density_log = em
    variable = DiffusiveFlux_em
    block = 0
    position_units = ${dom0Scale}
  []
  [EFieldAdvAux_emliq]
    type = ADEFieldAdvAux
    potential = potential
    density_log = emliq
    variable = EFieldAdvAux_emliq
    block = 1
    position_units = ${dom1Scale}
  []
  [DiffusiveFlux_emliq]
    type = ADDiffusiveFlux
    density_log = emliq
    variable = DiffusiveFlux_emliq
    block = 1
    position_units = ${dom1Scale}
  []
[]

[InterfaceKernels]
  [em_advection]
    type = InterfaceAdvection
    mean_en_neighbor = mean_en
    potential_neighbor = potential
    neighbor_var = em
    variable = emliq
    boundary = master1_interface
    position_units = ${dom1Scale}
    neighbor_position_units = ${dom0Scale}
  []
  [em_diffusion]
    type = InterfaceLogDiffusionElectrons
    mean_en_neighbor = mean_en
    neighbor_var = em
    variable = emliq
    boundary = master1_interface
    position_units = ${dom1Scale}
    neighbor_position_units = ${dom0Scale}
  []
[]

[BCs]
  [mean_en_physical_right]
    type = HagelaarEnergyBC
    variable = mean_en
    boundary = 'master0_interface'
    potential = potential
    electrons = em
    r = 0.99
    position_units = ${dom0Scale}
  []
  [mean_en_physical_left]
    type = HagelaarEnergyBC
    variable = mean_en
    boundary = 'left'
    potential = potential
    electrons = em
    r = 0
    position_units = ${dom0Scale}
  []
  [secondary_energy_left]
    type = SecondaryElectronEnergyBC
    variable = mean_en
    boundary = 'left'
    potential = potential
    electrons = em
    ions = 'Arp'
    r = 0
    emission_coeffs = 0.05
    position_units = ${dom0Scale}
    secondary_electron_energy = 3
  []

  [potential_left]
    type = NeumannCircuitVoltageMoles_KV
    variable = potential
    boundary = left
    function = potential_bc_func
    ions = Arp
    data_provider = data_provider
    electrons = em
    electron_energy = mean_en
    r = 0
    position_units = ${dom0Scale}
    emission_coeffs = 0.05
  []
  [potential_dirichlet_right]
    type = DirichletBC
    variable = potential
    boundary = right
    value = 0
  []
  [em_physical_right]
    type = HagelaarElectronBC
    variable = em
    boundary = 'master0_interface'
    potential = potential
    electron_energy = mean_en
    r = 0.99
    position_units = ${dom0Scale}
  []
  [Arp_physical_right_diffusion]
    type = HagelaarIonDiffusionBC
    variable = Arp
    boundary = 'master0_interface'
    r = 0
    position_units = ${dom0Scale}
  []
  [Arp_physical_right_advection]
    type = HagelaarIonAdvectionBC
    variable = Arp
    boundary = 'master0_interface'
    potential = potential
    r = 0
    position_units = ${dom0Scale}
  []

  [em_physical_left]
    type = HagelaarElectronBC
    variable = em
    boundary = 'left'
    potential = potential
    electron_energy = mean_en
    r = 0
    position_units = ${dom0Scale}
  []
  [sec_electrons_left]
    type = SecondaryElectronBC
    variable = em
    boundary = 'left'
    potential = potential
    ions = Arp
    electron_energy = mean_en
    r = 0
    position_units = ${dom0Scale}
    emission_coeffs = 0.05
  []
  [Arp_physical_left_diffusion]
    type = HagelaarIonDiffusionBC
    variable = Arp
    boundary = 'left'
    r = 0
    position_units = ${dom0Scale}
  []
  [Arp_physical_left_advection]
    type = HagelaarIonAdvectionBC
    variable = Arp
    boundary = 'left'
    potential = potential
    r = 0
    position_units = ${dom0Scale}
  []

  [emliq_right]
    type = DCIonBC
    variable = emliq
    boundary = right
    potential = potential
    position_units = ${dom1Scale}
  []
  [OHm_physical]
    type = DCIonBC
    variable = OHm
    boundary = 'right'
    potential = potential
    position_units = ${dom1Scale}
  []
[]

[ICs]
  [em_ic]
    type = ConstantIC
    variable = em
    value = -21
    block = 0
  []
  [emliq_ic]
    type = ConstantIC
    variable = emliq
    value = -21
    block = 1
  []
  [Arp_ic]
    type = ConstantIC
    variable = Arp
    value = -21
    block = 0
  []
  [mean_en_ic]
    type = ConstantIC
    variable = mean_en
    value = -20
    block = 0
  []
  [OHm_ic]
    type = ConstantIC
    variable = OHm
    value = -15.6
    block = 1
  []
  [potential_ic]
    type = FunctionIC
    variable = potential
    function = potential_ic_func
  []
[]

[Functions]
  [potential_bc_func]
    type = ParsedFunction
    # expression = '1.25*tanh(1e6*t)'
    expression = -1.25
  []
  [potential_ic_func]
    type = ParsedFunction
    expression = '-1.25 * (1.0001e-3 - x)'
  []
[]

[Materials]
  [gas_block]
    type = Gas
    interp_trans_coeffs = true
    interp_elastic_coeff = true
    ramp_trans_coeffs = false
    em = em
    potential = potential
    ip = Arp
    mean_en = mean_en
    user_se_coeff = 0.05
    block = 0
    property_tables_file = td_argon_mean_en.txt
  []
  [water_block]
    type = Water
    block = 1
    potential = potential
  []
[]

(test/tests/1d_dc/mean_en_multi.i)

dom0Scale = 1e-3
dom1Scale = 1e-7

[GlobalParams]
  offset = 20
  potential_units = kV
  use_moles = true
  # potential_units = V
[]

[Mesh]
  [file]
    type = FileMeshGenerator
    file = 'liquidNew.msh'
  []
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = '0'
    paired_block = '1'
    new_boundary = 'master0_interface'
    input = file
  []
  [interface_again]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = '1'
    paired_block = '0'
    new_boundary = 'master1_interface'
    input = interface
  []
  [left]
    type = SideSetsFromNormalsGenerator
    normals = '-1 0 0'
    new_boundary = 'left'
    input = interface_again
  []
  [right]
    type = SideSetsFromNormalsGenerator
    normals = '1 0 0'
    new_boundary = 'right'
    input = left
  []
[]

[Problem]
  type = FEProblem
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
    ksp_norm = none
  []
[]

[Executioner]
  type = Transient
  end_time = 1e-1
  petsc_options = '-snes_converged_reason -snes_linesearch_monitor'
  # petsc_options = '-snes_test_display'
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount -ksp_type -snes_linesearch_minlambda'
  petsc_options_value = 'lu NONZERO 1.e-10 preonly 1e-3'
  line_search = 'bt'
  nl_rel_tol = 1e-4
  nl_abs_tol = 7.6e-5
  dtmin = 1e-12
  [TimeSteppers]
    [Adaptive]
      type = IterationAdaptiveDT
      cutback_factor = 0.4
      dt = 1e-11
      growth_factor = 1.2
      optimal_iterations = 15
    []
  []
[]

[Outputs]
  perf_graph = true
  # print_linear_residuals = false
  [out]
    type = Exodus
    execute_on = 'final'
  []
  [dof_map]
    type = DOFMap
  []
[]

[Debug]
  show_var_residual_norms = true
[]

[UserObjects]
  [data_provider]
    type = ProvideMobility
    electrode_area = 5.02e-7 # Formerly 3.14e-6
    ballast_resist = 1e6
    e = 1.6e-19
  []
[]

[Kernels]
  [em_time_deriv]
    type = ElectronTimeDerivative
    variable = em
    block = 0
  []
  [em_advection]
    type = EFieldAdvection
    variable = em
    potential = potential
    block = 0
    position_units = ${dom0Scale}
  []
  [em_diffusion]
    type = CoeffDiffusion
    variable = em
    block = 0
    position_units = ${dom0Scale}
  []
  [em_ionization]
    type = ElectronsFromIonization
    variable = em
    potential = potential
    mean_en = mean_en
    em = em
    block = 0
    position_units = ${dom0Scale}
  []
  [em_log_stabilization]
    type = LogStabilizationMoles
    variable = em
    block = 0
  []
  [emliq_time_deriv]
    type = ElectronTimeDerivative
    variable = emliq
    block = 1
  []
  [emliq_advection]
    type = EFieldAdvection
    variable = emliq
    potential = potential
    block = 1
    position_units = ${dom1Scale}
  []
  [emliq_diffusion]
    type = CoeffDiffusion
    variable = emliq
    block = 1
    position_units = ${dom1Scale}
  []
  [emliq_reactant_first_order_rxn]
    type = ReactantFirstOrderRxn
    variable = emliq
    block = 1
  []
  [emliq_water_bi_sink]
    type = ReactantAARxn
    variable = emliq
    block = 1
  []
  [emliq_log_stabilization]
    type = LogStabilizationMoles
    variable = emliq
    block = 1
  []
  [potential_diffusion_dom1]
    type = CoeffDiffusionLin
    variable = potential
    block = 0
    position_units = ${dom0Scale}
  []
  [potential_diffusion_dom2]
    type = CoeffDiffusionLin
    variable = potential
    block = 1
    position_units = ${dom1Scale}
  []
  [Arp_charge_source]
    type = ChargeSourceMoles_KV
    variable = potential
    charged = Arp
    block = 0
  []
  [em_charge_source]
    type = ChargeSourceMoles_KV
    variable = potential
    charged = em
    block = 0
  []
  [emliq_charge_source]
    type = ChargeSourceMoles_KV
    variable = potential
    charged = emliq
    block = 1
  []
  [OHm_charge_source]
    type = ChargeSourceMoles_KV
    variable = potential
    charged = OHm
    block = 1
  []

  [Arp_time_deriv]
    type = ElectronTimeDerivative
    variable = Arp
    block = 0
  []

  [Arp_advection]
    type = EFieldAdvection
    variable = Arp
    potential = potential
    position_units = ${dom0Scale}
    block = 0
  []

  [ArEx_excitation]
    type = ExcitationReaction
    variable = ArEx
    potential = potential
    em = em
    mean_en = mean_en
    block = 0
    position_units = ${dom0Scale}
    reactant = false
  []

  [ArEx_deexcitation]
    type = ExcitationReaction
    variable = ArEx
    potential = potential
    em = em
    mean_en = mean_en
    block = 0
    position_units = ${dom0Scale}
    reactant = true
  []

  [Arp_diffusion]
    type = CoeffDiffusion
    variable = Arp
    block = 0
    position_units = ${dom0Scale}
  []
  [Arp_ionization]
    type = IonsFromIonization
    variable = Arp
    potential = potential
    em = em
    mean_en = mean_en
    block = 0
    position_units = ${dom0Scale}
  []
  [Arp_log_stabilization]
    type = LogStabilizationMoles
    variable = Arp
    block = 0
  []

  [ArEx_time_deriv]
    type = ElectronTimeDerivative
    variable = ArEx
    block = 0
  []

  [ArEx_diffusion]
    type = CoeffDiffusion
    variable = ArEx
    block = 0
    position_units = ${dom0Scale}
  []

  [ArEx_log_stabilization]
    type = LogStabilizationMoles
    variable = ArEx
    block = 0
    offset = 30
  []

  [ArTest_time_deriv]
    type = ElectronTimeDerivative
    variable = ArTest
    block = 0
  []

  [ArTest_diffusion]
    type = CoeffDiffusion
    variable = ArTest
    block = 0
    position_units = ${dom0Scale}
  []

  [ArTest_log_stabilization]
    type = LogStabilizationMoles
    variable = ArTest
    block = 0
  []

  [OHm_time_deriv]
    type = ElectronTimeDerivative
    variable = OHm
    block = 1
  []
  [OHm_advection]
    type = EFieldAdvection
    variable = OHm
    potential = potential
    block = 1
    position_units = ${dom1Scale}
  []
  [OHm_diffusion]
    type = CoeffDiffusion
    variable = OHm
    block = 1
    position_units = ${dom1Scale}
  []
  [OHm_log_stabilization]
    type = LogStabilizationMoles
    variable = OHm
    block = 1
  []
  [OHm_product_first_order_rxn]
    type = ProductFirstOrderRxn
    variable = OHm
    v = emliq
    block = 1
  []
  [OHm_product_aabb_rxn]
    type = ProductAABBRxn
    variable = OHm
    v = emliq
    block = 1
  []

  [mean_en_time_deriv]
    type = ElectronTimeDerivative
    variable = mean_en
    block = 0
  []
  [mean_en_advection]
    type = EFieldAdvection
    variable = mean_en
    potential = potential
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_diffusion]
    type = CoeffDiffusion
    variable = mean_en
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_joule_heating]
    type = JouleHeating
    variable = mean_en
    potential = potential
    em = em
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_ionization]
    type = ElectronEnergyLossFromIonization
    variable = mean_en
    potential = potential
    em = em
    block = 0
    position_units = ${dom0Scale}
  []

  [mean_en_elastic]
    type = ElectronEnergyLossFromElastic
    variable = mean_en
    potential = potential
    em = em
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_excitation]
    type = ElectronEnergyLossFromExcitation
    variable = mean_en
    potential = potential
    em = em
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_log_stabilization]
    type = LogStabilizationMoles
    variable = mean_en
    block = 0
    offset = 15
  []
[]

[Variables]
  [potential]
  []
  [em]
    block = 0
  []
  [emliq]
    block = 1
    # scaling = 1e-5
  []

  [Arp]
    block = 0
  []

  [ArEx]
    block = 0
  []

  [ArTest]
    block = 0
  []

  [mean_en]
    block = 0
    # scaling = 1e-1
  []

  [OHm]
    block = 1
    # scaling = 1e-5
  []
[]

[AuxVariables]
  [e_temp]
    block = 0
    order = CONSTANT
    family = MONOMIAL
  []
  [x]
    order = CONSTANT
    family = MONOMIAL
  []
  [x_node]
  []
  [rho]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [rholiq]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [em_lin]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [emliq_lin]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [Arp_lin]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [ArEx_lin]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [ArTest_lin]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [OHm_lin]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [Efield]
    order = CONSTANT
    family = MONOMIAL
  []
  [Current_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [Current_emliq]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [Current_Arp]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [Current_OHm]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [tot_gas_current]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [tot_liq_current]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [tot_flux_OHm]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [EFieldAdvAux_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [DiffusiveFlux_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [EFieldAdvAux_emliq]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [DiffusiveFlux_emliq]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [PowerDep_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [PowerDep_Arp]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [ProcRate_el]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [ProcRate_ex]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [ProcRate_iz]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
[]

[AuxKernels]
  [PowerDep_em]
    type = ADPowerDep
    density_log = em
    potential = potential
    art_diff = false
    potential_units = kV
    variable = PowerDep_em
    position_units = ${dom0Scale}
    block = 0
  []
  [PowerDep_Arp]
    type = ADPowerDep
    density_log = Arp
    potential = potential
    art_diff = false
    potential_units = kV
    variable = PowerDep_Arp
    position_units = ${dom0Scale}
    block = 0
  []
  [ProcRate_el]
    type = ADProcRate
    em = em
    potential = potential
    proc = el
    variable = ProcRate_el
    position_units = ${dom0Scale}
    block = 0
  []
  [ProcRate_ex]
    type = ADProcRate
    em = em
    potential = potential
    proc = ex
    variable = ProcRate_ex
    position_units = ${dom0Scale}
    block = 0
  []
  [ProcRate_iz]
    type = ADProcRate
    em = em
    potential = potential
    proc = iz
    variable = ProcRate_iz
    position_units = ${dom0Scale}
    block = 0
  []
  [e_temp]
    type = ElectronTemperature
    variable = e_temp
    electron_density = em
    mean_en = mean_en
    block = 0
  []
  [x_g]
    type = Position
    variable = x
    position_units = ${dom0Scale}
    block = 0
  []
  [x_l]
    type = Position
    variable = x
    position_units = ${dom1Scale}
    block = 1
  []
  [x_ng]
    type = Position
    variable = x_node
    position_units = ${dom0Scale}
    block = 0
  []
  [x_nl]
    type = Position
    variable = x_node
    position_units = ${dom1Scale}
    block = 1
  []
  [rho]
    type = ParsedAux
    variable = rho
    coupled_variables = 'em_lin Arp_lin'
    expression = 'Arp_lin - em_lin'
    execute_on = 'timestep_end'
    block = 0
  []
  [rholiq]
    type = ParsedAux
    variable = rholiq
    coupled_variables = 'emliq_lin OHm_lin' # H3Op_lin OHm_lin'
    expression = '-emliq_lin - OHm_lin' # 'H3Op_lin - em_lin - OHm_lin'
    execute_on = 'timestep_end'
    block = 1
  []
  [tot_gas_current]
    type = ParsedAux
    variable = tot_gas_current
    coupled_variables = 'Current_em Current_Arp'
    expression = 'Current_em + Current_Arp'
    execute_on = 'timestep_end'
    block = 0
  []
  [tot_liq_current]
    type = ParsedAux
    variable = tot_liq_current
    coupled_variables = 'Current_emliq Current_OHm' # Current_H3Op Current_OHm'
    expression = 'Current_emliq + Current_OHm' # + Current_H3Op + Current_OHm'
    execute_on = 'timestep_end'
    block = 1
  []
  [em_lin]
    type = DensityMoles
    variable = em_lin
    density_log = em
    block = 0
  []
  [emliq_lin]
    type = DensityMoles
    variable = emliq_lin
    density_log = emliq
    block = 1
  []
  [Arp_lin]
    type = DensityMoles
    variable = Arp_lin
    density_log = Arp
    block = 0
  []
  [ArEx_lin]
    type = DensityMoles
    variable = ArEx_lin
    density_log = ArEx
    block = 0
  []
  [ArTest_lin]
    type = DensityMoles
    variable = ArTest_lin
    density_log = ArTest
    block = 0
  []
  [OHm_lin]
    type = DensityMoles
    variable = OHm_lin
    density_log = OHm
    block = 1
  []
  [Efield_g]
    type = Efield
    component = 0
    potential = potential
    variable = Efield
    position_units = ${dom0Scale}
    block = 0
  []
  [Efield_l]
    type = Efield
    component = 0
    potential = potential
    variable = Efield
    position_units = ${dom1Scale}
    block = 1
  []
  [Current_em]
    type = ADCurrent
    potential = potential
    density_log = em
    variable = Current_em
    art_diff = false
    block = 0
    position_units = ${dom0Scale}
  []
  [Current_emliq]
    type = ADCurrent
    potential = potential
    density_log = emliq
    variable = Current_emliq
    art_diff = false
    block = 1
    position_units = ${dom1Scale}
  []
  [Current_Arp]
    type = ADCurrent
    potential = potential
    density_log = Arp
    variable = Current_Arp
    art_diff = false
    block = 0
    position_units = ${dom0Scale}
  []
  [Current_OHm]
    block = 1
    type = ADCurrent
    potential = potential
    density_log = OHm
    variable = Current_OHm
    art_diff = false
    position_units = ${dom1Scale}
  []
  [tot_flux_OHm]
    block = 1
    type = ADTotalFlux
    potential = potential
    density_log = OHm
    variable = tot_flux_OHm
  []
  [EFieldAdvAux_em]
    type = ADEFieldAdvAux
    potential = potential
    density_log = em
    variable = EFieldAdvAux_em
    block = 0
    position_units = ${dom0Scale}
  []
  [DiffusiveFlux_em]
    type = ADDiffusiveFlux
    density_log = em
    variable = DiffusiveFlux_em
    block = 0
    position_units = ${dom0Scale}
  []
  [EFieldAdvAux_emliq]
    type = ADEFieldAdvAux
    potential = potential
    density_log = emliq
    variable = EFieldAdvAux_emliq
    block = 1
    position_units = ${dom1Scale}
  []
  [DiffusiveFlux_emliq]
    type = ADDiffusiveFlux
    density_log = emliq
    variable = DiffusiveFlux_emliq
    block = 1
    position_units = ${dom1Scale}
  []
[]

[InterfaceKernels]
  [em_advection]
    type = InterfaceAdvection
    mean_en_neighbor = mean_en
    potential_neighbor = potential
    neighbor_var = em
    variable = emliq
    boundary = master1_interface
    position_units = ${dom1Scale}
    neighbor_position_units = ${dom0Scale}
  []
  [em_diffusion]
    type = InterfaceLogDiffusionElectrons
    mean_en_neighbor = mean_en
    neighbor_var = em
    variable = emliq
    boundary = master1_interface
    position_units = ${dom1Scale}
    neighbor_position_units = ${dom0Scale}
  []
[]

[BCs]
  [potential_left]
    type = NeumannCircuitVoltageMoles_KV
    variable = potential
    boundary = left
    function = potential_bc_func
    ions = Arp
    data_provider = data_provider
    electrons = em
    electron_energy = mean_en
    r = 0
    position_units = ${dom0Scale}
    emission_coeffs = 0.05
  []
  [potential_dirichlet_right]
    type = DirichletBC
    variable = potential
    boundary = right
    value = 0
  []
  [em_physical_right]
    type = HagelaarElectronBC
    variable = em
    boundary = 'master0_interface'
    potential = potential
    electron_energy = mean_en
    r = 0.99
    position_units = ${dom0Scale}
  []
  [Arp_physical_right_diffusion]
    type = HagelaarIonDiffusionBC
    variable = Arp
    boundary = 'master0_interface'
    r = 0
    position_units = ${dom0Scale}
  []
  [Arp_physical_right_advection]
    type = HagelaarIonAdvectionBC
    variable = Arp
    boundary = 'master0_interface'
    potential = potential
    r = 0
    position_units = ${dom0Scale}
  []
  [mean_en_physical_right]
    type = HagelaarEnergyBC
    variable = mean_en
    boundary = 'master0_interface'
    potential = potential
    electrons = em
    r = 0.99
    position_units = ${dom0Scale}
  []
  [em_physical_left]
    type = HagelaarElectronBC
    variable = em
    boundary = 'left'
    potential = potential
    electron_energy = mean_en
    r = 0
    position_units = ${dom0Scale}
  []
  [sec_electrons_left]
    type = SecondaryElectronBC
    variable = em
    boundary = 'left'
    potential = potential
    ions = Arp
    electron_energy = mean_en
    r = 0
    position_units = ${dom0Scale}
    emission_coeffs = 0.05
  []
  [Arp_physical_left_diffusion]
    type = HagelaarIonDiffusionBC
    variable = Arp
    boundary = 'left'
    r = 0
    position_units = ${dom0Scale}
  []
  [Arp_physical_left_advection]
    type = HagelaarIonAdvectionBC
    variable = Arp
    boundary = 'left'
    potential = potential
    r = 0
    position_units = ${dom0Scale}
  []
  [ArEx_physical_right_diffusion]
    type = HagelaarIonDiffusionBC
    variable = ArEx
    boundary = 'master0_interface'
    r = 0
    position_units = ${dom0Scale}
  []
  [ArEx_physical_left_diffusion]
    type = HagelaarIonDiffusionBC
    variable = ArEx
    boundary = 'left'
    r = 0
    position_units = ${dom0Scale}
  []
  [ArTest_physical_right_diffusion]
    type = HagelaarIonDiffusionBC
    variable = ArTest
    boundary = 'master0_interface'
    r = 0
    position_units = ${dom0Scale}
  []
  [ArTest_physical_left_diffusion]
    type = HagelaarIonDiffusionBC
    variable = ArTest
    boundary = 'left'
    r = 0
    position_units = ${dom0Scale}
  []
  [mean_en_physical_left]
    type = HagelaarEnergyBC
    variable = mean_en
    boundary = 'left'
    potential = potential
    electrons = em
    r = 0
    position_units = ${dom0Scale}
  []
  [secondary_energy_left]
    type = SecondaryElectronEnergyBC
    variable = mean_en
    boundary = 'left'
    potential = potential
    electrons = em
    ions = 'Arp'
    r = 0
    emission_coeffs = 0.05
    secondary_electron_energy = 3
    position_units = ${dom0Scale}
  []
  [emliq_right]
    type = DCIonBC
    variable = emliq
    boundary = right
    potential = potential
    position_units = ${dom1Scale}
  []
  [OHm_physical]
    type = DCIonBC
    variable = OHm
    boundary = 'right'
    potential = potential
    position_units = ${dom1Scale}
  []
[]

[ICs]
  [em_ic]
    type = ConstantIC
    variable = em
    value = -21
    block = 0
  []
  [emliq_ic]
    type = ConstantIC
    variable = emliq
    value = -21
    block = 1
  []
  [Arp_ic]
    type = ConstantIC
    variable = Arp
    value = -21
    block = 0
  []
  [ArEx_ic]
    type = ConstantIC
    variable = ArEx
    value = -16
    #value = -21
    block = 0
  []
  [ArTest_ic]
    type = ConstantIC
    variable = ArTest
    value = -31
    #value = -21
    block = 0
  []
  [mean_en_ic]
    type = ConstantIC
    variable = mean_en
    value = -20
    block = 0
  []
  [potential_ic]
    type = FunctionIC
    variable = potential
    function = potential_ic_func
  []
  [OHm_ic]
    type = ConstantIC
    variable = OHm
    value = -15.6
    block = 1
  []
[]

[Functions]
  [potential_bc_func]
    type = ParsedFunction
    # expression = '1.25*tanh(1e6*t)'
    expression = -1.25
  []
  [potential_ic_func]
    type = ParsedFunction
    expression = '-1.25 * (1.0001e-3 - x)'
  []
[]

[Materials]
  [gas_block]
    type = GasBase
    interp_trans_coeffs = true
    interp_elastic_coeff = true
    ramp_trans_coeffs = false
    em = em
    potential = potential
    ip = Arp
    mean_en = mean_en
    user_se_coeff = .05
    block = 0
    property_tables_file = td_argon_mean_en.txt
    position_units = ${dom0Scale}
  []
  [gas_species_0]
    type = ADHeavySpecies
    heavy_species_name = Arp
    heavy_species_mass = 6.64e-26
    heavy_species_charge = 1.0
    block = 0
  []
  [gas_species_1]
    type = ADHeavySpecies
    heavy_species_name = ArEx
    heavy_species_mass = 6.64e-26
    heavy_species_charge = 1.0
    block = 0
  []
  [gas_species_2]
    type = ADHeavySpecies
    heavy_species_name = ArTest
    heavy_species_mass = 6.64e-26
    heavy_species_charge = 1.0
    block = 0
  []
  [water_block]
    type = Water
    block = 1
    potential = potential
  []
[]

(test/tests/crane_action/townsend_units.i)

# THIS INPUT FILE IS BASED ON mean_en.i (1d_dc test)

dom0Scale = 1e-3
dom1Scale = 1e-7

[GlobalParams]
  offset = 20
  potential_units = kV
  use_moles = true
[]

[Mesh]
  [file]
    type = FileMeshGenerator
    file = 'townsend_units.msh'
  []
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = '0'
    paired_block = '1'
    new_boundary = 'master0_interface'
    input = file
  []
  [interface_again]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = '1'
    paired_block = '0'
    new_boundary = 'master1_interface'
    input = interface
  []
  [left]
    type = SideSetsFromNormalsGenerator
    normals = '-1 0 0'
    new_boundary = 'left'
    input = interface_again
  []
  [right]
    type = SideSetsFromNormalsGenerator
    normals = '1 0 0'
    new_boundary = 'right'
    input = left
  []
[]

[Problem]
  type = FEProblem
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  end_time = 1e-1
  petsc_options = '-snes_converged_reason -snes_linesearch_monitor'
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount'
  petsc_options_value = 'lu NONZERO 1.e-10'
  nl_rel_tol = 1e-4
  nl_abs_tol = 7.6e-5
  dtmin = 1e-15
  l_max_its = 20
  [TimeSteppers]
    [Adaptive]
      type = IterationAdaptiveDT
      cutback_factor = 0.4
      dt = 1e-11
      growth_factor = 1.2
      optimal_iterations = 30
    []
  []
[]

[Outputs]
  perf_graph = true
  [out]
    type = Exodus
    execute_on = 'final'
  []
  #[dof_map]
  #  type = DOFMap
  #[]
[]

[Debug]
  #show_var_residual_norms = true
[]

[UserObjects]
  [data_provider]
    type = ProvideMobility
    electrode_area = 5.02e-7 # Formerly 3.14e-6
    ballast_resist = 1e6
    e = 1.6e-19
  []
[]

[Kernels]
  [em_time_deriv]
    type = ElectronTimeDerivative
    variable = em
    block = 0
  []
  [em_advection]
    type = EFieldAdvection
    variable = em
    potential = potential
    block = 0
    position_units = ${dom0Scale}
  []
  [em_diffusion]
    type = CoeffDiffusion
    variable = em
    block = 0
    position_units = ${dom0Scale}
  []

  [em_log_stabilization]
    type = LogStabilizationMoles
    variable = em
    block = 0
  []

  [emliq_time_deriv]
    type = ElectronTimeDerivative
    variable = emliq
    block = 1
  []
  [emliq_advection]
    type = EFieldAdvection
    variable = emliq
    potential = potential
    block = 1
    position_units = ${dom1Scale}
  []
  [emliq_diffusion]
    type = CoeffDiffusion
    variable = emliq
    block = 1
    position_units = ${dom1Scale}
  []

  [emliq_log_stabilization]
    type = LogStabilizationMoles
    variable = emliq
    block = 1
  []

  [potential_diffusion_dom1]
    type = CoeffDiffusionLin
    variable = potential
    block = 0
    position_units = ${dom0Scale}
  []
  [potential_diffusion_dom2]
    type = CoeffDiffusionLin
    variable = potential
    block = 1
    position_units = ${dom1Scale}
  []
  [Arp_charge_source]
    type = ChargeSourceMoles_KV
    variable = potential
    charged = Arp
    block = 0
  []
  [em_charge_source]
    type = ChargeSourceMoles_KV
    variable = potential
    charged = em
    block = 0
  []
  [emliq_charge_source]
    type = ChargeSourceMoles_KV
    variable = potential
    charged = emliq
    block = 1
  []
  [OHm_charge_source]
    type = ChargeSourceMoles_KV
    variable = potential
    charged = OHm
    block = 1
  []

  [Arp_time_deriv]
    type = ElectronTimeDerivative
    variable = Arp
    block = 0
  []
  [Arp_advection]
    type = EFieldAdvection
    variable = Arp
    potential = potential
    position_units = ${dom0Scale}
    block = 0
  []
  [Arp_diffusion]
    type = CoeffDiffusion
    variable = Arp
    block = 0
    position_units = ${dom0Scale}
  []

  [Arp_log_stabilization]
    type = LogStabilizationMoles
    variable = Arp
    block = 0
  []

  [OHm_time_deriv]
    type = ElectronTimeDerivative
    variable = OHm
    block = 1
  []
  [OHm_advection]
    type = EFieldAdvection
    variable = OHm
    potential = potential
    block = 1
    position_units = ${dom1Scale}
  []
  [OHm_diffusion]
    type = CoeffDiffusion
    variable = OHm
    block = 1
    position_units = ${dom1Scale}
  []
  [OHm_log_stabilization]
    type = LogStabilizationMoles
    variable = OHm
    block = 1
  []

  [mean_en_time_deriv]
    type = ElectronTimeDerivative
    variable = mean_en
    block = 0
  []
  [mean_en_advection]
    type = EFieldAdvection
    variable = mean_en
    potential = potential
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_diffusion]
    type = CoeffDiffusion
    variable = mean_en
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_joule_heating]
    type = JouleHeating
    variable = mean_en
    potential = potential
    em = em
    block = 0
    position_units = ${dom0Scale}
  []
  [mean_en_log_stabilization]
    type = LogStabilizationMoles
    variable = mean_en
    block = 0
    offset = 15
  []
[]

[Variables]
  [potential]
  []

  [em]
    block = 0
  []

  [emliq]
    block = 1
  []

  [Arp]
    block = 0
  []

  [mean_en]
    block = 0
  []

  [OHm]
    block = 1
  []
[]

[AuxVariables]
  [H2O]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 10.92252
    block = 1
  []
  [Ar]
    block = 0
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 3.70109
  []
  [e_temp]
    block = 0
    order = CONSTANT
    family = MONOMIAL
  []
  [x]
    order = CONSTANT
    family = MONOMIAL
  []
  [x_node]
  []
  [rho]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [rholiq]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [em_lin]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [emliq_lin]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [Arp_lin]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [OHm_lin]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [Efield]
    order = CONSTANT
    family = MONOMIAL
  []
  [Current_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [Current_emliq]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [Current_Arp]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [Current_OHm]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [tot_gas_current]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [tot_liq_current]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [tot_flux_OHm]
    block = 1
    order = CONSTANT
    family = MONOMIAL
  []
  [EFieldAdvAux_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [DiffusiveFlux_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [EFieldAdvAux_emliq]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [DiffusiveFlux_emliq]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [PowerDep_em]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [PowerDep_Arp]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  #[ProcRate_el]
  # order = CONSTANT
  # family = MONOMIAL
  # block = 0
  #[]
  #[ProcRate_ex]
  # order = CONSTANT
  # family = MONOMIAL
  # block = 0
  #[]
  #[ProcRate_iz]
  # order = CONSTANT
  # family = MONOMIAL
  # block = 0
  #[]
[]

[AuxKernels]
  [PowerDep_em]
    type = ADPowerDep
    density_log = em
    potential = potential
    art_diff = false
    potential_units = kV
    variable = PowerDep_em
    position_units = ${dom0Scale}
    block = 0
  []
  [PowerDep_Arp]
    type = ADPowerDep
    density_log = Arp
    potential = potential
    art_diff = false
    potential_units = kV
    variable = PowerDep_Arp
    position_units = ${dom0Scale}
    block = 0
  []
  #[ProcRate_el]
  #  type = ProcRate
  #  em = em
  #  potential = potential
  #  proc = el
  #  variable = ProcRate_el
  #  position_units = ${dom0Scale}
  #  block = 0
  #[]
  #[ProcRate_ex]
  #  type = ProcRate
  #  em = em
  #  potential = potential
  #  proc = ex
  #  variable = ProcRate_ex
  #  position_units = ${dom0Scale}
  #  block = 0
  #[]
  #[ProcRate_iz]
  #  type = ProcRate
  #  em = em
  #  potential = potential
  #  proc = iz
  #  variable = ProcRate_iz
  #  position_units = ${dom0Scale}
  #  block = 0
  #[]
  [e_temp]
    type = ElectronTemperature
    variable = e_temp
    electron_density = em
    mean_en = mean_en
    block = 0
  []
  [x_g]
    type = Position
    variable = x
    position_units = ${dom0Scale}
    block = 0
  []
  [x_l]
    type = Position
    variable = x
    position_units = ${dom1Scale}
    block = 1
  []
  [x_ng]
    type = Position
    variable = x_node
    position_units = ${dom0Scale}
    block = 0
  []
  [x_nl]
    type = Position
    variable = x_node
    position_units = ${dom1Scale}
    block = 1
  []
  [rho]
    type = ParsedAux
    variable = rho
    coupled_variables = 'em_lin Arp_lin'
    expression = 'Arp_lin - em_lin'
    execute_on = 'timestep_end'
    block = 0
  []
  [rholiq]
    type = ParsedAux
    variable = rholiq
    coupled_variables = 'emliq_lin OHm_lin' # H3Op_lin OHm_lin'
    expression = '-emliq_lin - OHm_lin' # 'H3Op_lin - em_lin - OHm_lin'
    execute_on = 'timestep_end'
    block = 1
  []
  [tot_gas_current]
    type = ParsedAux
    variable = tot_gas_current
    coupled_variables = 'Current_em Current_Arp'
    expression = 'Current_em + Current_Arp'
    execute_on = 'timestep_end'
    block = 0
  []
  [tot_liq_current]
    type = ParsedAux
    variable = tot_liq_current
    coupled_variables = 'Current_emliq Current_OHm' # Current_H3Op Current_OHm'
    expression = 'Current_emliq + Current_OHm' # + Current_H3Op + Current_OHm'
    execute_on = 'timestep_end'
    block = 1
  []
  [em_lin]
    type = DensityMoles
    variable = em_lin
    density_log = em
    block = 0
  []
  [emliq_lin]
    type = DensityMoles
    variable = emliq_lin
    density_log = emliq
    block = 1
  []
  [Arp_lin]
    type = DensityMoles
    variable = Arp_lin
    density_log = Arp
    block = 0
  []
  [OHm_lin]
    type = DensityMoles
    variable = OHm_lin
    density_log = OHm
    block = 1
  []
  [Efield_g]
    type = Efield
    component = 0
    potential = potential
    variable = Efield
    position_units = ${dom0Scale}
    block = 0
  []
  [Efield_l]
    type = Efield
    component = 0
    potential = potential
    variable = Efield
    position_units = ${dom1Scale}
    block = 1
  []
  [Current_em]
    type = ADCurrent
    potential = potential
    density_log = em
    variable = Current_em
    art_diff = false
    block = 0
    position_units = ${dom0Scale}
  []
  [Current_emliq]
    type = ADCurrent
    potential = potential
    density_log = emliq
    variable = Current_emliq
    art_diff = false
    block = 1
    position_units = ${dom1Scale}
  []
  [Current_Arp]
    type = ADCurrent
    potential = potential
    density_log = Arp
    variable = Current_Arp
    art_diff = false
    block = 0
    position_units = ${dom0Scale}
  []
  [Current_OHm]
    block = 1
    type = ADCurrent
    potential = potential
    density_log = OHm
    variable = Current_OHm
    art_diff = false
    position_units = ${dom1Scale}
  []
  [tot_flux_OHm]
    block = 1
    type = ADTotalFlux
    potential = potential
    density_log = OHm
    variable = tot_flux_OHm
  []
  [EFieldAdvAux_em]
    type = ADEFieldAdvAux
    potential = potential
    density_log = em
    variable = EFieldAdvAux_em
    block = 0
    position_units = ${dom0Scale}
  []
  [DiffusiveFlux_em]
    type = ADDiffusiveFlux
    density_log = em
    variable = DiffusiveFlux_em
    block = 0
    position_units = ${dom0Scale}
  []
  [EFieldAdvAux_emliq]
    type = ADEFieldAdvAux
    potential = potential
    density_log = emliq
    variable = EFieldAdvAux_emliq
    block = 1
    position_units = ${dom1Scale}
  []
  [DiffusiveFlux_emliq]
    type = ADDiffusiveFlux
    density_log = emliq
    variable = DiffusiveFlux_emliq
    block = 1
    position_units = ${dom1Scale}
  []
[]

[InterfaceKernels]
  [em_advection]
    type = InterfaceAdvection
    mean_en_neighbor = mean_en
    potential_neighbor = potential
    neighbor_var = em
    variable = emliq
    boundary = master1_interface
    position_units = ${dom1Scale}
    neighbor_position_units = ${dom0Scale}
  []
  [em_diffusion]
    type = InterfaceLogDiffusionElectrons
    mean_en_neighbor = mean_en
    neighbor_var = em
    variable = emliq
    boundary = master1_interface
    position_units = ${dom1Scale}
    neighbor_position_units = ${dom0Scale}
  []
[]

[BCs]
  [mean_en_physical_right]
    type = HagelaarEnergyBC
    variable = mean_en
    boundary = 'master0_interface'
    potential = potential
    electrons = em
    r = 0.99
    #r = 0.0
    position_units = ${dom0Scale}
  []
  [mean_en_physical_left]
    type = HagelaarEnergyBC
    variable = mean_en
    boundary = 'left'
    potential = potential
    electrons = em
    r = 0
    position_units = ${dom0Scale}
  []
  [secondary_energy_left]
    type = SecondaryElectronEnergyBC
    variable = mean_en
    boundary = 'left'
    potential = potential
    electrons = em
    ions = 'Arp'
    r = 0
    emission_coeffs = 0.05
    secondary_electron_energy = 3
    position_units = ${dom0Scale}
  []

  [potential_left]
    type = NeumannCircuitVoltageMoles_KV
    variable = potential
    boundary = left
    function = potential_bc_func
    ions = Arp
    data_provider = data_provider
    electrons = em
    electron_energy = mean_en
    r = 0
    position_units = ${dom0Scale}
    emission_coeffs = 0.05
  []
  [potential_dirichlet_right]
    type = DirichletBC
    variable = potential
    boundary = right
    value = 0
  []
  [em_physical_right]
    type = HagelaarElectronBC
    variable = em
    boundary = 'master0_interface'
    potential = potential
    electron_energy = mean_en
    r = 0.99
    #r = 0.0
    position_units = ${dom0Scale}
  []
  [Arp_physical_right_diffusion]
    type = HagelaarIonDiffusionBC
    variable = Arp
    boundary = 'master0_interface'
    r = 0
    position_units = ${dom0Scale}
  []
  [Arp_physical_right_advection]
    type = HagelaarIonAdvectionBC
    variable = Arp
    boundary = 'master0_interface'
    potential = potential
    r = 0
    position_units = ${dom0Scale}
  []

  [em_physical_left]
    type = HagelaarElectronBC
    variable = em
    boundary = 'left'
    potential = potential
    electron_energy = mean_en
    r = 0
    position_units = ${dom0Scale}
  []
  [sec_electrons_left]
    type = SecondaryElectronBC
    variable = em
    boundary = 'left'
    potential = potential
    ions = Arp
    electron_energy = mean_en
    r = 0
    position_units = ${dom0Scale}
    emission_coeffs = 0.05
  []
  [Arp_physical_left_diffusion]
    type = HagelaarIonDiffusionBC
    variable = Arp
    boundary = 'left'
    r = 0
    position_units = ${dom0Scale}
  []
  [Arp_physical_left_advection]
    type = HagelaarIonAdvectionBC
    variable = Arp
    boundary = 'left'
    potential = potential
    r = 0
    position_units = ${dom0Scale}
  []

  [emliq_right]
    type = DCIonBC
    variable = emliq
    boundary = right
    potential = potential
    position_units = ${dom1Scale}
  []
  [OHm_physical]
    type = DCIonBC
    variable = OHm
    boundary = 'right'
    potential = potential
    position_units = ${dom1Scale}
  []
[]

[ICs]
  [em_ic]
    type = ConstantIC
    variable = em
    value = -21
    block = 0
  []
  [emliq_ic]
    type = ConstantIC
    variable = emliq
    value = -21
    block = 1
  []
  [Arp_ic]
    type = ConstantIC
    variable = Arp
    value = -21
    block = 0
  []
  [mean_en_ic]
    type = ConstantIC
    variable = mean_en
    value = -20
    block = 0
  []
  [OHm_ic]
    type = ConstantIC
    variable = OHm
    value = -15.6
    block = 1
  []
  [potential_ic]
    type = FunctionIC
    variable = potential
    function = potential_ic_func
  []
[]

[Functions]
  [potential_bc_func]
    type = ParsedFunction
    # expression = '1.25*tanh(1e6*t)'
    expression = -1.25
  []
  [potential_ic_func]
    type = ParsedFunction
    expression = '-1.25 * (1.0001e-3 - x)'
  []
[]

[Materials]
  [water_block]
    type = Water
    block = 1
    potential = potential
  []

  [test]
    type = GasElectronMoments
    interp_trans_coeffs = true
    interp_elastic_coeff = true
    ramp_trans_coeffs = false
    user_p_gas = 101325
    user_se_coeff = 0.05
    em = em
    potential = potential
    mean_en = mean_en
    block = 0
    property_tables_file = 'townsend_coefficients/moments.txt'
  []

  [test_block1]
    type = GenericConstantMaterial
    block = 1
    prop_names = 'T_gas p_gas'
    prop_values = '300 1.01e5'
  []

  [gas_species_0]
    type = ADHeavySpecies
    heavy_species_name = Arp
    heavy_species_mass = 6.64e-26
    heavy_species_charge = 1.0
    block = 0
  []

  [gas_species_2]
    type = ADHeavySpecies
    heavy_species_name = Ar
    heavy_species_mass = 6.64e-26
    heavy_species_charge = 0.0
    block = 0
  []
[]

[Reactions]
  [Argon]
    species = 'em Arp'
    aux_species = 'Ar'
    reaction_coefficient_format = 'townsend'
    gas_species = 'Ar'
    electron_energy = 'mean_en'
    electron_density = 'em'
    include_electrons = true
    file_location = 'townsend_coefficients'
    potential = 'potential'
    use_log = true
    use_ad = true
    position_units = ${dom0Scale}
    track_rates = false
    block = 0

    reactions = 'em + Ar -> em + Ar               : EEDF [elastic] (reaction1)
                 em + Ar -> em + Ar*              : EEDF [-11.5] (reaction2)
                 em + Ar -> em + em + Arp         : EEDF [-15.76] (reaction3)'
  []

  [Water]
    species = 'emliq OHm'
    reaction_coefficient_format = 'rate'
    use_log = true
    use_ad = true
    aux_species = 'H2O'
    block = 1
    reactions = 'emliq -> H + OHm : 1064
                 emliq + emliq -> H2 + OHm + OHm : 3.136e8'
  []
[]

Overview
Example Input File Syntax
Input Parameters
Input Files