PotentialSurfaceCharge

Enforces the dielectric boundary condition on a potential variable at an interface, where the surface charge is provided as an ADMaterialProperty.

Overview

PotentialSurfaceCharge enforces the dielectric boundary condition at an interface by equating the difference in the electric field at an interface to a surface charge. To supply the surface charge, the ADSurfaceCharge object is often used.

The dielectric boundary condition is defined as

$var element = document.getElementById("moose-equation-f7afc2d5-c4e1-4335-bd23-6c661d5b7321");katex.render("\\left( \\varepsilon_{1} \\vec{E}_{1} - \\varepsilon_{2} \\vec{E}_{2} \\right) \\cdot \\textbf{n} = \\sigma", element, {displayMode:true,throwOnError:false});$

Where:

Subscript $var element = document.getElementById("moose-equation-d7928b36-85cd-4573-b87b-c04280109210");katex.render("1", element, {displayMode:false,throwOnError:false});$ is the the primary side of the interface,
Subscript $var element = document.getElementById("moose-equation-9a5adba7-9b81-415f-b746-a72875b6f0c0");katex.render("2", element, {displayMode:false,throwOnError:false});$ is the the neighboring side of the interface,
$var element = document.getElementById("moose-equation-c1d00943-dffd-4efe-ac93-d8d27b6ed057");katex.render("\\textbf{n}", element, {displayMode:false,throwOnError:false});$ is the normal vector of the interface,
$var element = document.getElementById("moose-equation-7e1c146a-24c3-46c0-bb75-29d73c4d557a");katex.render("\\vec{E}", element, {displayMode:false,throwOnError:false});$ is the electric field, and
$var element = document.getElementById("moose-equation-420ca241-521a-4647-b89d-fdf7ff385949");katex.render("\\sigma", element, {displayMode:false,throwOnError:false});$ is the surface charge.

PotentialSurfaceCharge assumes the electrostatic approximation (i.e., $var element = document.getElementById("moose-equation-ab252e55-9b6a-44c4-890c-faf74f2ec5a1");katex.render("\\vec{E} = -\\nabla V", element, {displayMode:false,throwOnError:false});$ ). When converting to the electrostatic potential and applying a scaling factor of the mesh, PotentialSurfaceCharge is defined as

$var element = document.getElementById("moose-equation-05ef5dc9-e3fd-4473-9919-6e9818973f08");katex.render("\\left( \\varepsilon_{2} \\left( \\nabla V_{2} / l_{c,2} \\right) - \\varepsilon_{1} \\nabla \\left( V_{1} / l_{c,1} \\right) \\right) \\cdot \\textbf{n} = \\sigma", element, {displayMode:true,throwOnError:false});$

Where:

$var element = document.getElementById("moose-equation-ac5337e6-ec89-4dda-a7ed-6dd15be64a05");katex.render("V", element, {displayMode:false,throwOnError:false});$ is the electrostatic potential, and
$var element = document.getElementById("moose-equation-fc6a43e4-9e56-46c5-8f5e-ec046ba61556");katex.render("l_{c}", element, {displayMode:false,throwOnError:false});$ is the scaling factor of the mesh.

Example Input File Syntax

[InterfaceKernels<<<{"href": "../../syntax/InterfaceKernels/index.html"}>>>]
  [potential_right]
    type = PotentialSurfaceCharge<<<{"description": "Enforces the dielectric boundary condition on a potential variable at an interface, where the surface charge is provided as an ADMaterialProperty.", "href": "PotentialSurfaceCharge.html"}>>>
    neighbor_var<<<{"description": "The variable on the other side of the interface."}>>> = potential_dom1
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = potential_dom0
    position_units<<<{"description": "The units of position."}>>> = ${dom0Scale}
    neighbor_position_units<<<{"description": "The units of position."}>>> = ${dom1Scale}
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = plasma_right
  []
[]

Input Parameters

neighbor_varThe variable on the other side of the interface.
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The variable on the other side of the interface.
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

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
matrix_onlyFalseWhether this object is only doing assembly to matrices (no vectors)
Default:False
C++ Type:bool
Controllable:No
Description:Whether this object is only doing assembly to matrices (no vectors)
neighbor_position_units1The units of position.
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The units of position.
position_units1The units of position.
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The units of position.

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>
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>
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>
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
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
Options:nontime, time
Controllable:No
Description:The tag for the vectors this Kernel should fill

Contribution To Tagged Field Data 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
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

diag_save_inThe name of auxiliary variables to save this Kernel'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 Kernel's diagonal Jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
diag_save_in_var_sideThis parameter must exist if diag_save_in variables are specified and must have the same length as diag_save_in. This vector specifies whether the corresponding aux_var should save-in jacobian contributions from the primary ('p') or secondary side ('s').
C++ Type:MultiMooseEnum
Options:m, s
Controllable:No
Description:This parameter must exist if diag_save_in variables are specified and must have the same length as diag_save_in. This vector specifies whether the corresponding aux_var should save-in jacobian contributions from the primary ('p') or secondary side ('s').
save_inThe name of auxiliary variables to save this Kernel'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 Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
save_in_var_sideThis parameter must exist if save_in variables are specified and must have the same length as save_in. This vector specifies whether the corresponding aux_var should save-in residual contributions from the primary ('p') or secondary side ('s').
C++ Type:MultiMooseEnum
Options:m, s
Controllable:No
Description:This parameter must exist if save_in variables are specified and must have the same length as save_in. This vector specifies whether the corresponding aux_var should save-in residual contributions from the primary ('p') or secondary side ('s').

Residual And Jacobian Debug Output 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/interface_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 = ' e       N_A      k_boltz  eps         se_energy T_gas  massem   p_gas'
    prop_values = '1.6e-19 6.022e23 1.38e-23 8.854e-12   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 = ' e       N_A      k_boltz  eps         se_energy T_gas  massem   p_gas'
    prop_values = '1.6e-19 6.022e23 1.38e-23 8.854e-12   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/interface_test.i)

# This example is a simple steady state potential with a constant
# surface charge.
# The surface charge is taken into account with an interface kernel,
# causing a discontinuity in the electric field.
# A MatchedValueBC is also needed to ensure the potential is
# continuous across the two regions.

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

[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 = left dielectric
  # Block 1 = plasma region
  # Block 2 = right dielectric
  [file]
    type = FileMeshGenerator
    file = 'interface_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 = Steady
  automatic_scaling = true
  solve_type = NEWTON
  #petsc_options = '-snes_test_jacobian'
[]

[Outputs]
  perf_graph = true
  [out]
    type = Exodus
  []
[]

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

[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}
  []
[]

[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
  [potential_left]
    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
  []
[]

[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
  []
[]

[Functions]
  # Define a sinusoidal voltage pulse with a frequency of 50 kHz
  # Amplitude is set to 750 Volts
  # (Note that here the amplitude is set to 0.75. 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.75*sin(2*3.1415926*f0*t)'
    expression = '-0.75'
  []

  # 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.75 * (2.0001e-4 - x)'
  []
[]

[Materials]
  #########
  # Define secondary electron emission coefficients on the left and right
  # dielectrics.
  #########
  [se_left]
    type = GenericConstantMaterial
    boundary = 'plasma_right'
    prop_names = 'se_coeff'
    prop_values = '0.01'
  []

  [gas_constants]
    type = GenericConstantMaterial
    block = 1
    prop_names = ' e       N_A      k_boltz  eps         se_energy T_gas  massem   p_gas'
    prop_values = '1.6e-19 6.022e23 1.38e-23 8.854e-12   1.        400    9.11e-31 1.01e5'
  []

  # Create a constant surface charge on the boundary named 'plasma_side'
  [sc_test]
    type = ADGenericConstantMaterial
    boundary = 'plasma_right'
    prop_names = 'surface_charge'
    prop_values = '-1e-6'
  []

  [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-11'
    block = 1
  []

  [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