ADSurfaceCharge

Adds a surface charge material property based on the rate of change of the total charged flux to a boundary. (NOTE: this material is meant to be boundary-restricted.)

Overview

ADSurfaceCharge calculates the surface charge at a boundary. The surface charge material property is labeled surface_charge.

The surface charge is defined as

$var element = document.getElementById("moose-equation-4fe29bde-f21a-4da1-a6dc-04c64ec20b72");katex.render("\\sigma = \\int \\vec{J} \\ \\text{d}t \\\\[10pt] \\vec{J} \\cdot \\textbf{n} = \\sum{\\left( q_{j} \\vec{\\Gamma}_{j} \\cdot \\textbf{n} \\right)} \\\\[10pt] \\vec{\\Gamma}_{j} = \\text{sign}_{j} \\mu_{j} \\vec{E} n_{j} - D_{j} \\nabla (n_{j}) \\\\[10pt]", element, {displayMode:true,throwOnError:false});$

Where:

$var element = document.getElementById("moose-equation-f016c86a-4fbf-412f-951c-0033a02dbdfd");katex.render("\\sigma", element, {displayMode:false,throwOnError:false});$ is the surface charge,
$var element = document.getElementById("moose-equation-faba4f3e-3d28-43ac-9c67-fdd62f55c3bf");katex.render("\\vec{J}", element, {displayMode:false,throwOnError:false});$ is the total current density,
the subscript $var element = document.getElementById("moose-equation-60ad70e7-19de-48e4-aa8f-f00a2e7104e3");katex.render("j", element, {displayMode:false,throwOnError:false});$ represents properties of a plasma species,
$var element = document.getElementById("moose-equation-11910209-01e8-49f9-a1bc-a210057141da");katex.render("q_{j}", element, {displayMode:false,throwOnError:false});$ is the charge of the species,
$var element = document.getElementById("moose-equation-4a1ca7f5-fd70-4d1b-acd0-31f362c8ca2f");katex.render("\\vec{\\Gamma}", element, {displayMode:false,throwOnError:false});$ is the drift-diffusion flux,
$var element = document.getElementById("moose-equation-f9ef303d-a040-4535-afbb-fe56fbe4888a");katex.render("\\textbf{n}", element, {displayMode:false,throwOnError:false});$ is the outward pointing unit normal on the boundary,
$var element = document.getElementById("moose-equation-9c4f18f1-ac4c-43e9-8139-320844687e91");katex.render("\\text{sign}_{j}", element, {displayMode:false,throwOnError:false});$ indicates the advection behavior ( $var element = document.getElementById("moose-equation-badfafd1-f63b-4b40-86ad-bacad505942e");katex.render("\\text{+}1", element, {displayMode:false,throwOnError:false});$ for positively charged species and $var element = document.getElementById("moose-equation-73d15b33-43df-4f2a-9c16-744e0e35f7d6");katex.render("\\text{-}1", element, {displayMode:false,throwOnError:false});$ for negatively charged species),
$var element = document.getElementById("moose-equation-1e0f74cb-76fb-477b-931f-0dcf81990672");katex.render("\\mu_{j}", element, {displayMode:false,throwOnError:false});$ is the mobility coefficient,
$var element = document.getElementById("moose-equation-6a30d74c-5dd5-4d9a-821a-aca5e18e2199");katex.render("D_{j}", element, {displayMode:false,throwOnError:false});$ is the diffusion coefficient,
$var element = document.getElementById("moose-equation-5bb0663c-de40-4866-b441-d8d0df621c26");katex.render("\\vec{E}", element, {displayMode:false,throwOnError:false});$ is the electric field, and
$var element = document.getElementById("moose-equation-25b7452f-03ea-4a7f-b98a-de85afff6928");katex.render("n_{j}", element, {displayMode:false,throwOnError:false});$ is the species density.

Using the midpoint method for integration, the surface charge calculation becomes

$var element = document.getElementById("moose-equation-229d0236-0219-4e03-a9cb-63a5792372c3");katex.render("\\sigma_{t} = \\sigma_{t-1} + \\vec{J} \\cdot \\textbf{n} \\ \\text{d}t", element, {displayMode:true,throwOnError:false});$

Where:

$var element = document.getElementById("moose-equation-21f4949d-1c2d-427a-abdc-f4c6831e2dc2");katex.render("\\sigma_{t}", element, {displayMode:false,throwOnError:false});$ is the surface charge of the current time step,
$var element = document.getElementById("moose-equation-198274bd-7333-4f12-aece-a1176c599162");katex.render("\\sigma_{t-1}", element, {displayMode:false,throwOnError:false});$ is the surface of the previous time step, and
$var element = document.getElementById("moose-equation-791e41b8-78db-4267-9bdd-6a4b6b6cc5dc");katex.render("\\text{d}t", element, {displayMode:false,throwOnError:false});$ is the difference between time steps.

When converting the density to logarithmic form and applying a scaling factor of the mesh, the drift-diffusion flux of ADSurfaceCharge is redefined as

$var element = document.getElementById("moose-equation-f4e6413f-0cef-42b0-8259-733cced4bd9d");katex.render("\\vec{\\Gamma}_{j} = N_{A} \\left(\\text{sign}_{j} \\mu_{j} \\frac{\\vec{E}}{l_{c}} \\exp(N_{j}) - D_{j} \\exp(N_{j}) \\frac{\\nabla (N_{j})}{l_{c}} \\right)", element, {displayMode:true,throwOnError:false});$

Where:

$var element = document.getElementById("moose-equation-b21c4417-5c86-4a52-9699-3849c672aaa0");katex.render("N_{j}", element, {displayMode:false,throwOnError:false});$ is the molar density of the species in logarithmic form,
$var element = document.getElementById("moose-equation-f4f073d9-139a-4723-b9ab-a2dab4801a95");katex.render("N_{A}", element, {displayMode:false,throwOnError:false});$ is Avogadro's number, and
$var element = document.getElementById("moose-equation-be1814ed-ccb4-41dc-89a3-0f7c43738c72");katex.render("l_{c}", element, {displayMode:false,throwOnError:false});$ is the scaling factor of the mesh.

Example Input File Syntax

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  [surface_charge_material]
    type = ADSurfaceCharge<<<{"description": "Adds a surface charge material property based on the rate of change of the total charged flux to a boundary. (NOTE: this material is meant to be boundary-restricted.)", "href": "ADSurfaceCharge.html"}>>>
    species<<<{"description": "All of the charged species that interact with this boundary."}>>> = 'neg pos'
    position_units<<<{"description": "Units of position."}>>> = ${dom0Scale}
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = 'plasma_right'
  []
[]

Input Parameters

position_unitsUnits of position.
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Units of position.
potential_unitsThe potential units.
C++ Type:std::string
Controllable:No
Description:The potential units.
speciesAll of the charged species that interact with this boundary.
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:All of the charged species that interact with this boundary.

Required Parameters

bc_typeHagelaarThe name of the family of BCs being used in this model. Options: Hagelaar (DEFAULT), Sakiyama, Lymberopoulos.
Default:Hagelaar
C++ Type:std::string
Controllable:No
Description:The name of the family of BCs being used in this model. Options: Hagelaar (DEFAULT), Sakiyama, Lymberopoulos.
blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
boundaryThe list of boundaries (ids or names) from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundaries (ids or names) from the mesh where this object applies
computeTrueWhen false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
Default:True
C++ Type:bool
Controllable:No
Description:When false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
constant_onNONEWhen ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
Default:NONE
C++ Type:MooseEnum
Options:NONE, ELEMENT, SUBDOMAIN
Controllable:No
Description:When ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
declare_suffixAn optional suffix parameter that can be appended to any declared 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 declared properties. The suffix will be prepended with a '_' character.
field_property_namefield_solver_interface_propertyName of the solver interface material property.
Default:field_solver_interface_property
C++ Type:std::string
Controllable:No
Description:Name of the solver interface material property.
ksThe recombination coefficient (for Lymberopoulos-type BC)
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The recombination coefficient (for Lymberopoulos-type BC)
r_electron0The electron reflection coefficient on this boundary.
Default:0
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The electron reflection coefficient on this boundary.
r_ion0The ion reflection coefficient on this boundary.
Default:0
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The ion reflection coefficient on this boundary.
secondary_electronsTrueWhether or not to include secondary electron emission in the surface charge calculation. Note that this should be consistent with the selected boundary conditions; if a secondary electron BC is used on this boundary, this should be true. DEFAULT: true.
Default:True
C++ Type:bool
Controllable:No
Description:Whether or not to include secondary electron emission in the surface charge calculation. Note that this should be consistent with the selected boundary conditions; if a secondary electron BC is used on this boundary, this should be true. DEFAULT: true.

Optional Parameters

control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
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
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
search_methodnearest_node_connected_sidesChoice of search algorithm. All options begin by finding the nearest node in the primary boundary to a query point in the secondary boundary. In the default nearest_node_connected_sides algorithm, primary boundary elements are searched iff that nearest node is one of their nodes. This is fast to determine via a pregenerated node-to-elem map and is robust on conforming meshes. In the optional all_proximate_sides algorithm, primary boundary elements are searched iff they touch that nearest node, even if they are not topologically connected to it. This is more CPU-intensive but is necessary for robustness on any boundary surfaces which has disconnections (such as Flex IGA meshes) or non-conformity (such as hanging nodes in adaptively h-refined meshes).
Default:nearest_node_connected_sides
C++ Type:MooseEnum
Options:nearest_node_connected_sides, all_proximate_sides
Controllable:No
Description:Choice of search algorithm. All options begin by finding the nearest node in the primary boundary to a query point in the secondary boundary. In the default nearest_node_connected_sides algorithm, primary boundary elements are searched iff that nearest node is one of their nodes. This is fast to determine via a pregenerated node-to-elem map and is robust on conforming meshes. In the optional all_proximate_sides algorithm, primary boundary elements are searched iff they touch that nearest node, even if they are not topologically connected to it. This is more CPU-intensive but is necessary for robustness on any boundary surfaces which has disconnections (such as Flex IGA meshes) or non-conformity (such as hanging nodes in adaptively h-refined meshes).
seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
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
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

output_propertiesList of material properties, from this material, to output (outputs must also be defined to an output type)
C++ Type:std::vector<std::string>
Controllable:No
Description:List of material properties, from this material, to output (outputs must also be defined to an output type)
outputsnone Vector of output names where you would like to restrict the output of variables(s) associated with this object
Default:none
C++ Type:std::vector<OutputName>
Controllable:No
Description:Vector of output names where you would like to restrict the output of variables(s) associated with this object

Outputs Parameters

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.
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
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.

Material Property Retrieval Parameters

Input Files

(test/tests/surface_charge/dbd_test.i)

(test/tests/surface_charge/dbd_test.i)

# This example sets up a simple varying-voltage discharge with a
# dielectric at the right boundary.
#
# Two charged particle species are considered, creatively named 'pos'
# and 'neg'. They both have the same initial conditions and boundary
# conditions; the only difference is their charge.
# As the voltage varies, the particle flux to the right boundary
# (including both advection and diffusion) will switch polarity and
# the surface charge will become increasingly negative or positive.
#
# The left boundary is a driven electrode with a sinusoidal waveform.
# A dielectric boundary condition including surface charge accumulation
# is included on the right boundary.
#
# Note that without surface charge the electric field across both
# boundaries should be identical since their permittivities are
# both set to the same value.
# Surface charge shields the electric field in the gas region,
# preventing strong field buildup.

dom0Scale = 1e-4
dom1Scale = 1e-4

[GlobalParams]
  #offset = 40
  potential_units = kV
  use_moles = true
[]

[Mesh]
  # The mesh file generates the appropriate 1D mesh, but the interfaces
  # needed for the potential and surface charge have not been defined.
  # Here SideSetsBetweenSubdomainsGenerator is used to generate the
  # appropriate interfaces at the dielectrics.
  #
  # Block 0 = plasma region
  # Block 1 = dielectric
  [file]
    type = FileMeshGenerator
    file = 'dbd_mesh.msh'
  []
  [plasma_right]
    # plasma master
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = '0'
    paired_block = '1'
    new_boundary = 'plasma_right'
    input = file
  []
  [dielectric_left]
    # left dielectric master
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = '1'
    paired_block = '0'
    new_boundary = 'dielectric_left'
    input = plasma_right
  []
  [left]
    type = SideSetsFromNormalsGenerator
    normals = '-1 0 0'
    new_boundary = 'left'
    input = dielectric_left
  []
  [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 = 10e-5
  petsc_options = '-snes_converged_reason -snes_linesearch_monitor'
  #num_steps = 1
  #petsc_options = '-snes_converged_reason -snes_linesearch_monitor -snes_test_jacobian'
  solve_type = NEWTON
  line_search = 'basic'
  petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount -snes_stol'
  petsc_options_value = 'lu NONZERO 1.e-10 0'
  nl_rel_tol = 1e-8
  nl_div_tol = 1e4
  l_max_its = 100
  nl_max_its = 25

  dtmin = 1e-22
  dtmax = 1e-6
  #dt = 1e-8
  [TimeSteppers]
    [Adaptive]
      type = IterationAdaptiveDT
      cutback_factor = 0.4
      dt = 1e-9
      growth_factor = 1.2
      optimal_iterations = 30
    []
  []
[]

[Outputs]
  perf_graph = true
  [out]
    type = Exodus
    output_material_properties = true
    show_material_properties = 'surface_charge'
  []
[]

[AuxVariables]
  [neg_density]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [pos_density]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []

  [Efield]
    order = CONSTANT
    family = MONOMIAL
  []

  [x]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [neg_calc]
    type = DensityMoles
    variable = neg_density
    density_log = neg
    execute_on = 'initial timestep_end'
    block = 0
  []
  [pos_calc]
    type = DensityMoles
    variable = pos_density
    density_log = pos
    execute_on = 'initial timestep_end'
    block = 0
  []
  #[Arp_calc]
  #  type = DensityMoles
  #  variable = Arp_density
  #  density_log = Arp
  #  execute_on = 'initial timestep_end'
  #  block = 0
  #[]
  [x_d0]
    type = Position
    variable = x
    position_units = ${dom0Scale}
    execute_on = 'initial timestep_end'
    block = 0
  []
  [x_d1]
    type = Position
    variable = x
    position_units = ${dom1Scale}
    execute_on = 'initial timestep_end'
    block = 1
  []
  [Efield_g]
    type = Efield
    component = 0
    variable = Efield
    position_units = ${dom0Scale}
    block = 0
  []
  [Efield_l]
    type = Efield
    component = 0
    variable = Efield
    position_units = ${dom1Scale}
    block = 1
  []
[]

[Variables]
  [potential_dom0]
    block = 0
  []
  [potential_dom1]
    block = 1
  []

  [neg]
    block = 0
  []
  [pos]
    block = 0
  []
  #[mean_en]
  #  block = 0
  #[]
  #[Arp]
  #  block = 0
  #[]
[]

[Kernels]
  [potential_diffusion_dom0]
    type = CoeffDiffusionLin
    variable = potential_dom0
    block = 0
    position_units = ${dom0Scale}
  []
  [potential_diffusion_dom1]
    type = CoeffDiffusionLin
    variable = potential_dom1
    block = 1
    position_units = ${dom1Scale}
  []

  # Potential source terms
  #[em_source]
  #  type = ChargeSourceMoles_KV
  #  variable = potential_dom0
  #  charged = em
  #  block = 0
  #[]
  #[Arp_source]
  #  type = ChargeSourceMoles_KV
  #  variable = potential_dom0
  #  charged = Arp
  #  block = 0
  #[]
  #
  # Electrons
  #[d_em_dt]
  #  type = TimeDerivativeLog
  #  variable = em
  #  block = 0
  #[]
  [neg_advection]
    type = EFieldAdvection
    variable = neg
    position_units = ${dom0Scale}
    block = 0
  []
  [neg_diffusion]
    type = CoeffDiffusion
    variable = neg
    position_units = ${dom0Scale}
    block = 0
  []

  [pos_advection]
    type = EFieldAdvection
    variable = pos
    position_units = ${dom0Scale}
    block = 0
  []
  [pos_diffusion]
    type = CoeffDiffusion
    variable = pos
    position_units = ${dom0Scale}
    block = 0
  []
  #[em_offset]
  #  type = LogStabilizationMoles
  #  variable = em
  #  block = 0
  #[]
  #
  #[d_mean_en_dt]
  #  type = TimeDerivativeLog
  #  variable = mean_en
  #  block = 0
  #[]
  #[mean_en_advection]
  #  type = ADEFieldAdvection
  #  em = em
  #  variable = mean_en
  #  potential = potential_dom0
  #  position_units = ${dom0Scale}
  #  block = 0
  #[]
  #[mean_en_diffusion]
  #  type = ADCoeffDiffusion
  #  variable = mean_en
  #  position_units = ${dom0Scale}
  #  block = 0
  #[]
  #[mean_en_joule_heating]
  #  type = ADJouleHeating
  #  variable = mean_en
  #  em = em
  #  potential = potential_dom0
  #  position_units = ${dom0Scale}
  #  block = 0
  #[]
  #[mean_en_offset]
  #  type = LogStabilizationMoles
  #  variable = mean_en
  #  block = 0
  #[]
  #
  #[d_Arp_dt]
  #  type = TimeDerivativeLog
  #  variable = Arp
  #  block = 0
  #[]
  #[Arp_advection]
  #  type = ADEFieldAdvection
  #  variable = Arp
  #  potential = potential_dom0
  #  position_units = ${dom0Scale}
  #  block = 0
  #[]
  #[Arp_diffusion]
  #  type = ADCoeffDiffusion
  #  variable = Arp
  #  position_units = ${dom0Scale}
  #  block = 0
  #[]
  #[Arp_offset]
  #  type = LogStabilizationMoles
  #  variable = Arp
  #  block = 0
  #[]
[]

[InterfaceKernels]
  # At the dielectric interfaces, the potential is required to be continuous
  # in value but discontinuous in slope due to surface charge accumulation.
  #
  # The potential requires two different boundary conditions on each side:
  #    (1) An InterfaceKernel to provide the Neumann boundary condition
  #    (2) A MatchedValueBC to ensure that the potential remains continuous
  #
  # This interface kernel simply applies a dielectric BC with surface charge
  # accumulation.
  [potential_right]
    type = PotentialSurfaceCharge
    neighbor_var = potential_dom1
    variable = potential_dom0
    position_units = ${dom0Scale}
    neighbor_position_units = ${dom1Scale}
    boundary = plasma_right
  []
[]

[BCs]
  # Interface BCs:
  [match_potential]
    type = MatchedValueBC
    variable = potential_dom1
    v = potential_dom0
    boundary = 'dielectric_left'
  []

  # Electrode and ground BCs:
  # Electrode is at x = 0, named 'left'
  [potential_left]
    type = FunctionDirichletBC
    variable = potential_dom0
    function = potential_input
    boundary = 'left'
  []

  # Ground is at x = 0.3 mm, named 'right'
  [potential_dirichlet_right]
    type = DirichletBC
    variable = potential_dom1
    boundary = right
    value = 0
  []

  # Both charged species will have a specified dirichlet BC on the driven electrode
  # and a diffusion BC on the right side.
  # As the voltage varies with time, the charged particle flux on the right will
  # switch between positive and negative, effectively changing the surface
  # charge polarity.
  [neg_left]
    type = DirichletBC
    variable = neg
    value = -10
    boundary = 'left'
  []
  [neg_right]
    type = HagelaarIonDiffusionBC
    variable = neg
    r = 0
    position_units = ${dom0Scale}
    boundary = 'plasma_right'
  []

  [pos_left]
    type = DirichletBC
    variable = pos
    value = -10
    boundary = 'left'
  []
  [pos_right]
    type = HagelaarIonDiffusionBC
    variable = pos
    r = 0
    position_units = ${dom0Scale}
    boundary = 'plasma_right'
  []
[]

[ICs]
  [potential_dom0_ic]
    type = FunctionIC
    variable = potential_dom0
    function = potential_ic_func
  []
  [potential_dom1_ic]
    type = FunctionIC
    variable = potential_dom1
    function = potential_ic_func
  []

  [neg_ic]
    type = ConstantIC
    variable = neg
    value = -10
    block = 0
  []
  [pos_ic]
    type = ConstantIC
    variable = pos
    value = -10
    block = 0
  []
[]

[Functions]
  # Define a sinusoidal voltage pulse with a frequency of 50 kHz
  # Amplitude is set to 200 Volts
  # (Note that here the amplitude is set to 0.2. Potential units are
  # typically included as kV, not V. This option is set in the
  # GlobalParams block.)
  [potential_input]
    type = ParsedFunction
    symbol_names = 'f0'
    symbol_values = '50e3'
    expression = '-0.2*sin(2*3.1415926*f0*t)'
  []

  # Set the initial condition to a line from -10 V on the left and
  # 0 on the right.
  # (Poisson solver tends to struggle with a uniformly zero potential IC.)
  [potential_ic_func]
    type = ParsedFunction
    expression = '-0.001 * (2.000e-4 - x)'
  []
[]

[Materials]
  # This material creates a boundary-restricted material property called "surface_charge"
  # The value of surface charge is based on the total charged particle flux impacting
  # the specified boundary.
  # In this case we have two charged species, 'pos' (positively charged) and 'neg'
  # (negatively charged).
  [surface_charge_material]
    type = ADSurfaceCharge
    species = 'neg pos'
    position_units = ${dom0Scale}
    boundary = 'plasma_right'
  []

  # This just defines some constants that Zapdos needs to run.
  [gas_constants]
    type = GenericConstantMaterial
    block = 0
    prop_names = ' se_energy T_gas  massem   p_gas'
    prop_values = '1.        400    9.11e-31 1.01e5'
  []

  #
  [dielectric_left_side]
    type = ADGenericConstantMaterial
    prop_names = 'diffpotential_dom0'
    prop_values = '8.85e-12'
    block = 0
  []
  [gas_phase]
    type = ADGenericConstantMaterial
    prop_names = 'diffpotential_dom1'
    prop_values = '8.85e-12'
    block = 1
  []

  ######
  # HeavySpeciesMaterial defines charge, mass, transport coefficients, and
  # temperature for each species.
  #
  # Transport coefficients and temperature are defined as ADMaterialProperties.
  # Although they currently (as of June 16, 2020) remain constant, future
  # implementations may include mixture-averaged diffusion coefficients and
  # effective ion temperatures with nonlinear dependence on other variables.
  ######
  [gas_species_0]
    type = ADHeavySpecies
    heavy_species_name = neg
    heavy_species_mass = 6.64e-26
    heavy_species_charge = -1.0
    diffusivity = 1.6897e-5
    block = 0
  []
  [gas_species_2]
    type = ADHeavySpecies
    heavy_species_name = pos
    heavy_species_mass = 6.64e-26
    heavy_species_charge = 1.0
    diffusivity = 1.6897e-5
    block = 0
  []

  [field_solver0]
    type = FieldSolverMaterial
    potential = potential_dom0
    block = 0
  []

  [field_solver1]
    type = FieldSolverMaterial
    potential = potential_dom1
    block = 1
  []
[]

(test/tests/surface_charge/dbd_test.i)

# This example sets up a simple varying-voltage discharge with a
# dielectric at the right boundary.
#
# Two charged particle species are considered, creatively named 'pos'
# and 'neg'. They both have the same initial conditions and boundary
# conditions; the only difference is their charge.
# As the voltage varies, the particle flux to the right boundary
# (including both advection and diffusion) will switch polarity and
# the surface charge will become increasingly negative or positive.
#
# The left boundary is a driven electrode with a sinusoidal waveform.
# A dielectric boundary condition including surface charge accumulation
# is included on the right boundary.
#
# Note that without surface charge the electric field across both
# boundaries should be identical since their permittivities are
# both set to the same value.
# Surface charge shields the electric field in the gas region,
# preventing strong field buildup.

dom0Scale = 1e-4
dom1Scale = 1e-4

[GlobalParams]
  #offset = 40
  potential_units = kV
  use_moles = true
[]

[Mesh]
  # The mesh file generates the appropriate 1D mesh, but the interfaces
  # needed for the potential and surface charge have not been defined.
  # Here SideSetsBetweenSubdomainsGenerator is used to generate the
  # appropriate interfaces at the dielectrics.
  #
  # Block 0 = plasma region
  # Block 1 = dielectric
  [file]
    type = FileMeshGenerator
    file = 'dbd_mesh.msh'
  []
  [plasma_right]
    # plasma master
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = '0'
    paired_block = '1'
    new_boundary = 'plasma_right'
    input = file
  []
  [dielectric_left]
    # left dielectric master
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = '1'
    paired_block = '0'
    new_boundary = 'dielectric_left'
    input = plasma_right
  []
  [left]
    type = SideSetsFromNormalsGenerator
    normals = '-1 0 0'
    new_boundary = 'left'
    input = dielectric_left
  []
  [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 = 10e-5
  petsc_options = '-snes_converged_reason -snes_linesearch_monitor'
  #num_steps = 1
  #petsc_options = '-snes_converged_reason -snes_linesearch_monitor -snes_test_jacobian'
  solve_type = NEWTON
  line_search = 'basic'
  petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount -snes_stol'
  petsc_options_value = 'lu NONZERO 1.e-10 0'
  nl_rel_tol = 1e-8
  nl_div_tol = 1e4
  l_max_its = 100
  nl_max_its = 25

  dtmin = 1e-22
  dtmax = 1e-6
  #dt = 1e-8
  [TimeSteppers]
    [Adaptive]
      type = IterationAdaptiveDT
      cutback_factor = 0.4
      dt = 1e-9
      growth_factor = 1.2
      optimal_iterations = 30
    []
  []
[]

[Outputs]
  perf_graph = true
  [out]
    type = Exodus
    output_material_properties = true
    show_material_properties = 'surface_charge'
  []
[]

[AuxVariables]
  [neg_density]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []
  [pos_density]
    order = CONSTANT
    family = MONOMIAL
    block = 0
  []

  [Efield]
    order = CONSTANT
    family = MONOMIAL
  []

  [x]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [neg_calc]
    type = DensityMoles
    variable = neg_density
    density_log = neg
    execute_on = 'initial timestep_end'
    block = 0
  []
  [pos_calc]
    type = DensityMoles
    variable = pos_density
    density_log = pos
    execute_on = 'initial timestep_end'
    block = 0
  []
  #[Arp_calc]
  #  type = DensityMoles
  #  variable = Arp_density
  #  density_log = Arp
  #  execute_on = 'initial timestep_end'
  #  block = 0
  #[]
  [x_d0]
    type = Position
    variable = x
    position_units = ${dom0Scale}
    execute_on = 'initial timestep_end'
    block = 0
  []
  [x_d1]
    type = Position
    variable = x
    position_units = ${dom1Scale}
    execute_on = 'initial timestep_end'
    block = 1
  []
  [Efield_g]
    type = Efield
    component = 0
    variable = Efield
    position_units = ${dom0Scale}
    block = 0
  []
  [Efield_l]
    type = Efield
    component = 0
    variable = Efield
    position_units = ${dom1Scale}
    block = 1
  []
[]

[Variables]
  [potential_dom0]
    block = 0
  []
  [potential_dom1]
    block = 1
  []

  [neg]
    block = 0
  []
  [pos]
    block = 0
  []
  #[mean_en]
  #  block = 0
  #[]
  #[Arp]
  #  block = 0
  #[]
[]

[Kernels]
  [potential_diffusion_dom0]
    type = CoeffDiffusionLin
    variable = potential_dom0
    block = 0
    position_units = ${dom0Scale}
  []
  [potential_diffusion_dom1]
    type = CoeffDiffusionLin
    variable = potential_dom1
    block = 1
    position_units = ${dom1Scale}
  []

  # Potential source terms
  #[em_source]
  #  type = ChargeSourceMoles_KV
  #  variable = potential_dom0
  #  charged = em
  #  block = 0
  #[]
  #[Arp_source]
  #  type = ChargeSourceMoles_KV
  #  variable = potential_dom0
  #  charged = Arp
  #  block = 0
  #[]
  #
  # Electrons
  #[d_em_dt]
  #  type = TimeDerivativeLog
  #  variable = em
  #  block = 0
  #[]
  [neg_advection]
    type = EFieldAdvection
    variable = neg
    position_units = ${dom0Scale}
    block = 0
  []
  [neg_diffusion]
    type = CoeffDiffusion
    variable = neg
    position_units = ${dom0Scale}
    block = 0
  []

  [pos_advection]
    type = EFieldAdvection
    variable = pos
    position_units = ${dom0Scale}
    block = 0
  []
  [pos_diffusion]
    type = CoeffDiffusion
    variable = pos
    position_units = ${dom0Scale}
    block = 0
  []
  #[em_offset]
  #  type = LogStabilizationMoles
  #  variable = em
  #  block = 0
  #[]
  #
  #[d_mean_en_dt]
  #  type = TimeDerivativeLog
  #  variable = mean_en
  #  block = 0
  #[]
  #[mean_en_advection]
  #  type = ADEFieldAdvection
  #  em = em
  #  variable = mean_en
  #  potential = potential_dom0
  #  position_units = ${dom0Scale}
  #  block = 0
  #[]
  #[mean_en_diffusion]
  #  type = ADCoeffDiffusion
  #  variable = mean_en
  #  position_units = ${dom0Scale}
  #  block = 0
  #[]
  #[mean_en_joule_heating]
  #  type = ADJouleHeating
  #  variable = mean_en
  #  em = em
  #  potential = potential_dom0
  #  position_units = ${dom0Scale}
  #  block = 0
  #[]
  #[mean_en_offset]
  #  type = LogStabilizationMoles
  #  variable = mean_en
  #  block = 0
  #[]
  #
  #[d_Arp_dt]
  #  type = TimeDerivativeLog
  #  variable = Arp
  #  block = 0
  #[]
  #[Arp_advection]
  #  type = ADEFieldAdvection
  #  variable = Arp
  #  potential = potential_dom0
  #  position_units = ${dom0Scale}
  #  block = 0
  #[]
  #[Arp_diffusion]
  #  type = ADCoeffDiffusion
  #  variable = Arp
  #  position_units = ${dom0Scale}
  #  block = 0
  #[]
  #[Arp_offset]
  #  type = LogStabilizationMoles
  #  variable = Arp
  #  block = 0
  #[]
[]

[InterfaceKernels]
  # At the dielectric interfaces, the potential is required to be continuous
  # in value but discontinuous in slope due to surface charge accumulation.
  #
  # The potential requires two different boundary conditions on each side:
  #    (1) An InterfaceKernel to provide the Neumann boundary condition
  #    (2) A MatchedValueBC to ensure that the potential remains continuous
  #
  # This interface kernel simply applies a dielectric BC with surface charge
  # accumulation.
  [potential_right]
    type = PotentialSurfaceCharge
    neighbor_var = potential_dom1
    variable = potential_dom0
    position_units = ${dom0Scale}
    neighbor_position_units = ${dom1Scale}
    boundary = plasma_right
  []
[]

[BCs]
  # Interface BCs:
  [match_potential]
    type = MatchedValueBC
    variable = potential_dom1
    v = potential_dom0
    boundary = 'dielectric_left'
  []

  # Electrode and ground BCs:
  # Electrode is at x = 0, named 'left'
  [potential_left]
    type = FunctionDirichletBC
    variable = potential_dom0
    function = potential_input
    boundary = 'left'
  []

  # Ground is at x = 0.3 mm, named 'right'
  [potential_dirichlet_right]
    type = DirichletBC
    variable = potential_dom1
    boundary = right
    value = 0
  []

  # Both charged species will have a specified dirichlet BC on the driven electrode
  # and a diffusion BC on the right side.
  # As the voltage varies with time, the charged particle flux on the right will
  # switch between positive and negative, effectively changing the surface
  # charge polarity.
  [neg_left]
    type = DirichletBC
    variable = neg
    value = -10
    boundary = 'left'
  []
  [neg_right]
    type = HagelaarIonDiffusionBC
    variable = neg
    r = 0
    position_units = ${dom0Scale}
    boundary = 'plasma_right'
  []

  [pos_left]
    type = DirichletBC
    variable = pos
    value = -10
    boundary = 'left'
  []
  [pos_right]
    type = HagelaarIonDiffusionBC
    variable = pos
    r = 0
    position_units = ${dom0Scale}
    boundary = 'plasma_right'
  []
[]

[ICs]
  [potential_dom0_ic]
    type = FunctionIC
    variable = potential_dom0
    function = potential_ic_func
  []
  [potential_dom1_ic]
    type = FunctionIC
    variable = potential_dom1
    function = potential_ic_func
  []

  [neg_ic]
    type = ConstantIC
    variable = neg
    value = -10
    block = 0
  []
  [pos_ic]
    type = ConstantIC
    variable = pos
    value = -10
    block = 0
  []
[]

[Functions]
  # Define a sinusoidal voltage pulse with a frequency of 50 kHz
  # Amplitude is set to 200 Volts
  # (Note that here the amplitude is set to 0.2. Potential units are
  # typically included as kV, not V. This option is set in the
  # GlobalParams block.)
  [potential_input]
    type = ParsedFunction
    symbol_names = 'f0'
    symbol_values = '50e3'
    expression = '-0.2*sin(2*3.1415926*f0*t)'
  []

  # Set the initial condition to a line from -10 V on the left and
  # 0 on the right.
  # (Poisson solver tends to struggle with a uniformly zero potential IC.)
  [potential_ic_func]
    type = ParsedFunction
    expression = '-0.001 * (2.000e-4 - x)'
  []
[]

[Materials]
  # This material creates a boundary-restricted material property called "surface_charge"
  # The value of surface charge is based on the total charged particle flux impacting
  # the specified boundary.
  # In this case we have two charged species, 'pos' (positively charged) and 'neg'
  # (negatively charged).
  [surface_charge_material]
    type = ADSurfaceCharge
    species = 'neg pos'
    position_units = ${dom0Scale}
    boundary = 'plasma_right'
  []

  # This just defines some constants that Zapdos needs to run.
  [gas_constants]
    type = GenericConstantMaterial
    block = 0
    prop_names = ' se_energy T_gas  massem   p_gas'
    prop_values = '1.        400    9.11e-31 1.01e5'
  []

  #
  [dielectric_left_side]
    type = ADGenericConstantMaterial
    prop_names = 'diffpotential_dom0'
    prop_values = '8.85e-12'
    block = 0
  []
  [gas_phase]
    type = ADGenericConstantMaterial
    prop_names = 'diffpotential_dom1'
    prop_values = '8.85e-12'
    block = 1
  []

  ######
  # HeavySpeciesMaterial defines charge, mass, transport coefficients, and
  # temperature for each species.
  #
  # Transport coefficients and temperature are defined as ADMaterialProperties.
  # Although they currently (as of June 16, 2020) remain constant, future
  # implementations may include mixture-averaged diffusion coefficients and
  # effective ion temperatures with nonlinear dependence on other variables.
  ######
  [gas_species_0]
    type = ADHeavySpecies
    heavy_species_name = neg
    heavy_species_mass = 6.64e-26
    heavy_species_charge = -1.0
    diffusivity = 1.6897e-5
    block = 0
  []
  [gas_species_2]
    type = ADHeavySpecies
    heavy_species_name = pos
    heavy_species_mass = 6.64e-26
    heavy_species_charge = 1.0
    diffusivity = 1.6897e-5
    block = 0
  []

  [field_solver0]
    type = FieldSolverMaterial
    potential = potential_dom0
    block = 0
  []

  [field_solver1]
    type = FieldSolverMaterial
    potential = potential_dom1
    block = 1
  []
[]

Overview
Example Input File Syntax
Input Parameters
Input Files