DielectricBCWithEffEfield

Time integrated dielectric boundary condition (Based on Lymberopoulos and Economou (1994))

Overview

DielectricBCWithEffEfield is a type of NeumannBC for the potential on the boundary of a grounded ideal dielectric. The advection term for the ion is determined by an effective electric field.

The potential at the boundary of a grounded ideal dielectric is defined as

$var element = document.getElementById("moose-equation-78278a1c-c618-4730-a801-c9ddebd08fe1");katex.render("\\epsilon_{0}\\frac{\\partial (E \\cdot \\textbf{n}) }{\\partial t} - \\frac{\\epsilon_{i}}{d_{i}}\\frac{\\partial V_{i}}{\\partial t} = - e \\left( \\Gamma_{+} \\cdot \\textbf{n} -\\Gamma_{e} \\cdot \\textbf{n} \\right) \\\\[10pt] \\Gamma_{e} \\cdot \\textbf{n} = \\frac{1}{4}\\sqrt{\\frac{8 k T_{e}}{\\pi m_{e}}} \\ n_e - \\gamma \\Gamma_{+} \\cdot \\textbf{n} \\\\[10pt] \\Gamma_{+} \\cdot \\textbf{n} = a \\ \\mu_{+} \\ \\vec{E}_{\\text{Eff.}} \\cdot \\textbf{n} \\ n_{+} \\\\[10pt] a = \\begin{cases} 1, & \\mu_{+} \\ \\vec{E}_{\\text{Eff.}} \\cdot \\textbf{n} > 0\\\\ 0, & \\mu_{+} \\ \\vec{E}_{\\text{Eff.}} \\cdot \\textbf{n} \\leq 0\\\\ \\end{cases}", element, {displayMode:true,throwOnError:false});$

where

$var element = document.getElementById("moose-equation-e84e5357-0ea0-4767-8d65-a8a250ea47c6");katex.render("\\epsilon_{i}", element, {displayMode:false,throwOnError:false});$ is the permittivity of the dielectric,
$var element = document.getElementById("moose-equation-4dd2d8bf-4ee8-44cc-b449-b41ea4d7c1f0");katex.render("d_{i}", element, {displayMode:false,throwOnError:false});$ is the thickness of the dielectric,
$var element = document.getElementById("moose-equation-407838d1-1db1-41d5-91ec-711c0bbb7fba");katex.render("V_{i}", element, {displayMode:false,throwOnError:false});$ is the voltage on the dielectric,
$var element = document.getElementById("moose-equation-eb8e186c-a3ae-4164-911e-8cb6a3f517b4");katex.render("\\textbf{n}", element, {displayMode:false,throwOnError:false});$ is the normal to the boundary,
$var element = document.getElementById("moose-equation-1c89ed69-98ed-43c7-9eca-4ac49e679f0a");katex.render("e", element, {displayMode:false,throwOnError:false});$ is the elemental charge,
$var element = document.getElementById("moose-equation-f530f377-f14d-470e-bd62-3f3af9027dfc");katex.render("\\epsilon_{0}", element, {displayMode:false,throwOnError:false});$ is the permittivity of free space,
$var element = document.getElementById("moose-equation-0a819625-5f98-4196-85a1-7d8b8ca6fd9b");katex.render("\\vec{E} \\cdot \\textbf{n}", element, {displayMode:false,throwOnError:false});$ is the electric field normal to the dielectric,
$var element = document.getElementById("moose-equation-f85c7063-380e-49cd-a6e0-24180d9d1b7a");katex.render("\\Gamma_{e}", element, {displayMode:false,throwOnError:false});$ is the electron outflow flux, and
$var element = document.getElementById("moose-equation-86a27bd0-9e5c-490d-bf8d-964cfff87710");katex.render("\\Gamma_{+}", element, {displayMode:false,throwOnError:false});$ is the ion outflow flux.

note:Flux Information

$var element = document.getElementById("moose-equation-18d674b4-c8f0-47e6-b21e-39707175eff7");katex.render("\\Gamma_{e}", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-cea8a1dd-cab1-4a51-b48f-e81910919d21");katex.render("\\Gamma_{+}", element, {displayMode:false,throwOnError:false});$ are defined with the SakiyamaElectronDiffusionBC, SakiyamaSecondaryElectronWithEffEfieldBC and SakiyamaIonAdvectionWithEffEfieldBC (please refer to those BC's for more information on the fluxes).

To convert the above equation into the form of a NeumannBC, the time integral is taken such that:

$var element = document.getElementById("moose-equation-070bb146-edc5-4d0d-bddf-43bf7ca85af0");katex.render("\\int{ \\epsilon_{0}\\frac{\\partial (E \\cdot \\textbf{n}) }{\\partial t} } dt - \\int{ \\frac{\\epsilon_{i}}{d_{i}}\\frac{\\partial V_{i}}{\\partial t} } dt = \\int{ - e \\left( \\Gamma_{+} \\cdot \\textbf{n} -\\Gamma_{e} \\cdot \\textbf{n} \\right) } dt", element, {displayMode:true,throwOnError:false});$

Using the trapezoidal rule for the definite integral and rearranging the equation such that the electric field term is on one side results in:

$var element = document.getElementById("moose-equation-c2365753-38de-4678-975a-952b6b9247bb");katex.render("\\epsilon_{0} \\left( E \\cdot \\textbf{n} \\right) = \\epsilon_{0} \\left( E_{\\text{old}} \\cdot \\textbf{n} \\right) + \\frac{\\epsilon_{i}}{d_{i}} \\left( V - V_ {old} \\right) - 0.5 \\left( e \\left( \\Gamma_{+} \\cdot \\textbf{n} - \\Gamma_{e} \\cdot \\textbf{n} \\right) + e \\left( \\Gamma_{+} \\cdot \\textbf{n} - \\Gamma_{e} \\cdot \\textbf{n} \\right)_{old} \\right) dt", element, {displayMode:true,throwOnError:false});$

where

The subscript $var element = document.getElementById("moose-equation-68dbca5a-da46-479f-b614-5ea1fc05fcbe");katex.render("\\text{old}", element, {displayMode:false,throwOnError:false});$ is the value of the variable during the previous time step, and
$var element = document.getElementById("moose-equation-43431d8d-a072-47e3-b73d-70974f654211");katex.render("dt", element, {displayMode:false,throwOnError:false});$ is the current time step size.

Lastly, the electrostatic approximation is applied to the electric field normal to the dielectric, which results in a NeumannBC for the potential defined as:

$var element = document.getElementById("moose-equation-afabaf9c-0d3f-47df-bd0e-1aa6601f3c47");katex.render("\\epsilon_{0} \\left( -\\nabla V \\cdot \\textbf{n} \\right) = \\epsilon_{0} \\left( -\\nabla V_{\\text{old}} \\cdot \\textbf{n} \\right) + \\frac{\\epsilon_{i}}{d_{i}} \\left( V - V_ {old} \\right) - 0.5 \\left( e \\left( \\Gamma_{+} \\cdot \\textbf{n} - \\Gamma_{e} \\cdot \\textbf{n} \\right) + e \\left( \\Gamma_{+} \\cdot \\textbf{n} - \\Gamma_{e} \\cdot \\textbf{n} \\right)_{old} \\right) dt", element, {displayMode:true,throwOnError:false});$

Example Input File Syntax

[BCs]
  [potential_left_BC]
    type = DielectricBCWithEffEfield
    variable = potential
    electrons = em
    electron_energy = mean_en
    ions = 'ion'
    Ex = Ex
    Ey = Ey
    dielectric_constant = 1.0
    thickness = 1.0
    emission_coeffs = 'users_gamma'
    potential_units = V
    position_units = 1.0
    boundary = 3
  []
[]

Input Parameters

ExThe EField in the x-direction
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The EField in the x-direction
boundaryThe list of boundary IDs from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundary IDs from the mesh where this object applies
dielectric_constantThe dielectric constant of the material.
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The dielectric constant of the material.
electron_energyThe mean energy.
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The mean energy.
electronsThe electron density.
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The electron density.
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.
potential_unitsThe potential units.
C++ Type:std::string
Controllable:No
Description:The potential units.
thicknessThe thickness of the material.
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The thickness of the material.
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

EyThe EField in the y-direction
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The EField in the y-direction
EzThe EField in the z-direction
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The EField in the z-direction
displacementsThe displacements
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The displacements
emission_coeffsThe secondary electron emission coefficient for each ion provided in `ions`
C++ Type:std::vector<std::string>
Controllable:No
Description:The secondary electron emission coefficient for each ion provided in `ions`

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.
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
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
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
Controllable:No
Description:The seed for the master random number generator
skip_execution_outside_variable_domainFalseWhether to skip execution of this boundary condition when the variable it applies to is not defined on the boundary. This can facilitate setups with moving variable domains and fixed boundaries. Note that the FEProblem boundary-restricted integrity checks will also need to be turned off if using this option
Default:False
C++ Type:bool
Controllable:No
Description:Whether to skip execution of this boundary condition when the variable it applies to is not defined on the boundary. This can facilitate setups with moving variable domains and fixed boundaries. Note that the FEProblem boundary-restricted integrity checks will also need to be turned off if using this option
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

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/mms/bcs/2D_DielectricBCWithEffEfield.i)

References

Dimitris P Lymberopoulos and Demetre J Economou. Modeling and simulation of glow discharge plasma reactors. Journal of Vacuum Science & Technology A: Vacuum, Surfaces, and Films, 12(4):1229–1236, 1994. doi:10.1116/1.579300.

@article{lymberopoulos1994modeling,
    author = "Lymberopoulos, Dimitris P and Economou, Demetre J",
    title = "Modeling and simulation of glow discharge plasma reactors",
    journal = "Journal of Vacuum Science \\& Technology A: Vacuum, Surfaces, and Films",
    volume = "12",
    number = "4",
    pages = "1229--1236",
    year = "1994",
    publisher = "American Vacuum Society",
    doi = "10.1116/1.579300"
}

(test/tests/mms/bcs/2D_DielectricBCWithEffEfield.i)

[Mesh]
  [geo]
    type = FileMeshGenerator
    file = '2D_DielectricBCWithEffEfield_IC_out.e'
    use_for_exodus_restart = true
  []
[]

[Problem]
  type = FEProblem
[]

[Variables]
  [em]
    initial_from_file_var = em
  []
  [ion]
    initial_from_file_var = ion
  []
  [mean_en]
    initial_from_file_var = mean_en
  []

  [Ex]
    initial_from_file_var = Ex
  []
  [Ey]
    initial_from_file_var = Ey
  []

  [potential]
    initial_from_file_var = potential
  []
[]

# [ICs]
#   [em_IC]
#     type = FunctionIC
#     variable = em
#     function = 'em_ICs'
#   []
#   [ion_IC]
#     type = FunctionIC
#     variable = ion
#     function = 'ion_ICs'
#   []
#   [mean_en_IC]
#     type = FunctionIC
#     variable = mean_en
#     function = 'mean_en_ICs'
#   []
# []

[Kernels]
  #Electron Equations
  [em_time_derivative]
    type = ElectronTimeDerivative
    variable = em
  []
  [em_diffusion]
    type = CoeffDiffusion
    variable = em
    position_units = 1.0
  []
  [em_advection]
    type = EFieldAdvection
    variable = em
    position_units = 1.0
  []
  [em_source]
    type = BodyForce
    variable = em
    function = 'em_source'
  []

  #Ion Equations
  [ion_time_derivative]
    type = ElectronTimeDerivative
    variable = ion
  []
  [ion_diffusion]
    type = CoeffDiffusion
    variable = ion
    position_units = 1.0
  []
  [ion_advection]
    type = EffectiveEFieldAdvection
    variable = ion
    u = Ex
    v = Ey
    position_units = 1.0
  []
  [ion_source]
    type = BodyForce
    variable = ion
    function = 'ion_source'
  []

  #Eff. Efield
  [EffEfield_X_time_deriv]
    type = ADTimeDerivative
    variable = Ex
  []
  [EffEfield_X_diffusion]
    type = MatDiffusion
    D_name = diffEx
    variable = Ex
  []
  [EffEfield_X_source]
    type = BodyForce
    variable = Ex
    function = 'Ex_source'
  []
  [EffEfield_Y_time_deriv]
    type = ADTimeDerivative
    variable = Ey
  []
  [EffEfield_Y_diffusion]
    type = MatDiffusion
    D_name = diffEy
    variable = Ey
  []
  [EffEfield_Y_source]
    type = BodyForce
    variable = Ey
    function = 'Ey_source'
  []

  #Potential
  [Potential_time_deriv]
    type = ADTimeDerivative
    variable = potential
  []
  [Potential_diffusion]
    type = CoeffDiffusionLin
    variable = potential
    position_units = 1.0
  []
  [Potential_source]
    type = BodyForce
    variable = potential
    function = 'potential_source'
  []

  #Electron Energy Equations
  [mean_en_time_deriv]
    type = ElectronTimeDerivative
    variable = mean_en
  []
  [mean_en_diffusion]
    type = CoeffDiffusion
    variable = mean_en
    position_units = 1.0
  []
  [mean_en_source]
    type = BodyForce
    variable = mean_en
    function = 'energy_source'
  []
[]

[AuxVariables]
  [mean_en_sol]
  []

  [em_sol]
  []

  [ion_sol]
  []

  [Ex_sol]
  []
  [Ey_sol]
  []

  [potential_sol]
  []
[]

[AuxKernels]
  [mean_en_sol]
    type = FunctionAux
    variable = mean_en_sol
    function = mean_en_fun
  []

  [em_sol]
    type = FunctionAux
    variable = em_sol
    function = em_fun
  []

  [ion_sol]
    type = FunctionAux
    variable = ion_sol
    function = ion_fun
  []

  [Ex_sol]
    type = FunctionAux
    variable = Ex_sol
    function = Ex_fun
  []
  [Ey_sol]
    type = FunctionAux
    variable = Ey_sol
    function = Ey_fun
  []

  [potential_sol]
    type = FunctionAux
    variable = potential_sol
    function = potential_fun
  []
[]

[Functions]
  #Material Variables
  [massem]
    type = ConstantFunction
    value = 2.0
  []
  #Electron diffusion coeff.
  [diffem]
    type = ConstantFunction
    value = 0.05
  []
  [muem]
    type = ConstantFunction
    value = 0.01
  []
  #Electron energy mobility coeff.
  [diffmean_en]
    type = ConstantFunction
    value = 0.05
  []
  #Ion diffusion coeff.
  [diffion]
    type = ParsedFunction
    vars = diffem
    vals = diffem
    value = diffem
  []
  [muion]
    type = ParsedFunction
    vars = muem
    vals = muem
    value = muem
  []
  [N_A]
    type = ConstantFunction
    value = 1.0
  []
  [ee]
    type = ConstantFunction
    value = 1.0
  []
  [diffpotential]
    type = ConstantFunction
    value = 0.05
  []

  #Manufactured Solutions
  #The manufactured electron density solution
  [em_fun]
    type = ParsedFunction
    vars = 'N_A'
    vals = 'N_A'
    value = 'log((sin(pi*y) + 0.2*sin(2*pi*t)*cos(pi*y) + 2.0 + sin(pi*x)) / N_A)'
  []
  #The manufactured ion density solution
  [ion_fun]
    type = ParsedFunction
    vars = 'N_A'
    vals = 'N_A'
    value = 'log((sin(pi*y) + 0.2*sin(2*pi*t)*cos(pi*y) + 2.0 + sin(pi*x)) / N_A)'
  []
  #The manufactured electron energy solution
  [mean_en_fun]
    type = ParsedFunction
    vars = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    vals = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    value = 'log(((3*massem*pi*((diffpotential*pi*cos(pi*t) - cos(pi*t)*(sin(pi*x) + sin(pi*y)) + 4)/ee + 8*muion*pi*(sin(pi*t) + 1)*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2))^2)/(16*ee*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2))) / N_A)'
  []
  #The manufactured eff. Efield solution
  [Ex_fun]
    type = ParsedFunction
    vars = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    vals = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    value = '-pi*cos(pi*x)*(sin(pi*t) + 1)'
  []
  [Ey_fun]
    type = ParsedFunction
    vars = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    vals = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    value = '-pi*cos(pi*y)*(sin(pi*t) + 1)'
  []
  #The manufactured potential solution
  [potential_fun]
    type = ParsedFunction
    vars = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    vals = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    value = '0.25*(sin(pi*t)/pi + 1.0)*(sin(pi*y) + sin(pi*x)) - t'
  []

  #Source Terms in moles
  #The electron source term.
  [em_source]
    type = ParsedFunction
    vars = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    vals = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    value = '(diffem*pi^2*sin(pi*x) + (diffem*pi^2*(5*sin(pi*y) + cos(pi*y)*sin(2*pi*t)))/5 + (2*pi*cos(2*pi*t)*cos(pi*y))/5 - (muem*pi*(pi + sin(pi*t))*(10*sin(pi*x) + 10*sin(pi*y) + 10*sin(pi*x)*sin(pi*y) - 10*cos(pi*x)^2 - 10*cos(pi*y)^2 + cos(pi*y)*sin(2*pi*t)*sin(pi*x) + 2*cos(pi*y)*sin(2*pi*t)*sin(pi*y) + 10))/20) / N_A'
  []
  #The ion source term.
  [ion_source]
    type = ParsedFunction
    vars = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    vals = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    value = '(diffion*pi^2*sin(pi*x) + (diffion*pi^2*(5*sin(pi*y) + cos(pi*y)*sin(2*pi*t)))/5 + (2*pi*cos(2*pi*t)*cos(pi*y))/5 + (muion*pi^2*(sin(pi*t) + 1)*(10*sin(pi*x) + 10*sin(pi*y) + 10*sin(pi*x)*sin(pi*y) - 10*cos(pi*x)^2 - 10*cos(pi*y)^2 + cos(pi*y)*sin(2*pi*t)*sin(pi*x) + 2*cos(pi*y)*sin(2*pi*t)*sin(pi*y) + 10))/5) / N_A'
  []
  [energy_source]
    type = ParsedFunction
    vars = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    vals = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    value = '((3*massem*pi*((diffpotential*pi*cos(pi*t) - cos(pi*t)*(sin(pi*x) + sin(pi*y)) + 4)/ee + 8*muion*pi*(sin(pi*t) + 1)*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2))*(8*muion*pi^2*cos(pi*t)*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2) - (diffpotential*pi^2*sin(pi*t) - pi*sin(pi*t)*(sin(pi*x) + sin(pi*y)))/ee + (16*muion*pi^2*cos(2*pi*t)*cos(pi*y)*(sin(pi*t) + 1))/5))/(8*ee*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2)) - diffmean_en*((3*massem*pi*((pi*cos(pi*t)*cos(pi*y))/ee - 8*muion*pi*(pi*cos(pi*y) - (pi*sin(2*pi*t)*sin(pi*y))/5)*(sin(pi*t) + 1))^2)/(8*ee*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2)) + (3*massem*pi*((diffpotential*pi*cos(pi*t) - cos(pi*t)*(sin(pi*x) + sin(pi*y)) + 4)/ee + 8*muion*pi*(sin(pi*t) + 1)*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2))^2*(pi^2*sin(pi*y) + (pi^2*cos(pi*y)*sin(2*pi*t))/5))/(16*ee*(sin(x*pi) + sin(y*pi) + (cos(y*pi)*sin(2*t*pi))/5 + 2)^2) - (3*massem*pi*((diffpotential*pi*cos(pi*t) - cos(pi*t)*(sin(pi*x) + sin(pi*y)) + 4)/ee + 8*muion*pi*(sin(pi*t) + 1)*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2))*(8*muion*pi*(pi^2*sin(pi*y) + (pi^2*cos(pi*y)*sin(2*pi*t))/5)*(sin(pi*t) + 1) - (pi^2*cos(pi*t)*sin(pi*y))/ee))/(8*ee*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2)) + (3*massem*pi*((diffpotential*pi*cos(pi*t) - cos(pi*t)*(sin(pi*x) + sin(pi*y)) + 4)/ee + 8*muion*pi*(sin(pi*t) + 1)*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2))^2*(pi*cos(pi*y) - (pi*sin(2*pi*t)*sin(pi*y))/5)^2)/(8*ee*(sin(x*pi) + sin(y*pi) + (cos(y*pi)*sin(2*t*pi))/5 + 2)^3) + (3*massem*pi*((diffpotential*pi*cos(pi*t) - cos(pi*t)*(sin(pi*x) + sin(pi*y)) + 4)/ee + 8*muion*pi*(sin(pi*t) + 1)*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2))*((pi*cos(pi*t)*cos(pi*y))/ee - 8*muion*pi*(pi*cos(pi*y) - (pi*sin(2*pi*t)*sin(pi*y))/5)*(sin(pi*t) + 1))*(pi*cos(pi*y) - (pi*sin(2*pi*t)*sin(pi*y))/5))/(4*ee*(sin(x*pi) + sin(y*pi) + (cos(y*pi)*sin(2*t*pi))/5 + 2)^2)) - diffmean_en*((3*massem*pi*((pi*cos(pi*t)*cos(pi*x))/ee - 8*muion*pi^2*cos(pi*x)*(sin(pi*t) + 1))^2)/(8*ee*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2)) + (3*massem*pi^3*cos(pi*x)^2*((diffpotential*pi*cos(pi*t) - cos(pi*t)*(sin(pi*x) + sin(pi*y)) + 4)/ee + 8*muion*pi*(sin(pi*t) + 1)*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2))^2)/(8*ee*(sin(x*pi) + sin(y*pi) + (cos(y*pi)*sin(2*t*pi))/5 + 2)^3) + (3*massem*pi^3*sin(pi*x)*((diffpotential*pi*cos(pi*t) - cos(pi*t)*(sin(pi*x) + sin(pi*y)) + 4)/ee + 8*muion*pi*(sin(pi*t) + 1)*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2))^2)/(16*ee*(sin(x*pi) + sin(y*pi) + (cos(y*pi)*sin(2*t*pi))/5 + 2)^2) - (3*massem*pi*((diffpotential*pi*cos(pi*t) - cos(pi*t)*(sin(pi*x) + sin(pi*y)) + 4)/ee + 8*muion*pi*(sin(pi*t) + 1)*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2))*(8*muion*pi^3*sin(pi*x)*(sin(pi*t) + 1) - (pi^2*cos(pi*t)*sin(pi*x))/ee))/(8*ee*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2)) + (3*massem*pi^2*cos(pi*x)*((diffpotential*pi*cos(pi*t) - cos(pi*t)*(sin(pi*x) + sin(pi*y)) + 4)/ee + 8*muion*pi*(sin(pi*t) + 1)*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2))*((pi*cos(pi*t)*cos(pi*x))/ee - 8*muion*pi^2*cos(pi*x)*(sin(pi*t) + 1)))/(4*ee*(sin(x*pi) + sin(y*pi) + (cos(y*pi)*sin(2*t*pi))/5 + 2)^2)) - (3*massem*pi^2*cos(2*pi*t)*cos(pi*y)*((diffpotential*pi*cos(pi*t) - cos(pi*t)*(sin(pi*x) + sin(pi*y)) + 4)/ee + 8*muion*pi*(sin(pi*t) + 1)*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2))^2)/(40*ee*(sin(x*pi) + sin(y*pi) + (cos(y*pi)*sin(2*t*pi))/5 + 2)^2)) / N_A'
  []

  #The Ex source term.
  [Ex_source]
    type = ParsedFunction
    vars = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    vals = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    value = '-pi^2*cos(pi*t)*cos(pi*x) - diffpotential*pi^3*cos(pi*x)*(sin(pi*t) + 1)'
  []
  [Ey_source]
    type = ParsedFunction
    vars = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    vals = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    value = '-pi^2*cos(pi*t)*cos(pi*y) - diffpotential*pi^3*cos(pi*y)*(sin(pi*t) + 1)'
  []

  [potential_source]
    type = ParsedFunction
    vars = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    vals = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    value = '(cos(pi*t)*(sin(pi*x) + sin(pi*y)))/4 + diffpotential*pi^2*sin(pi*x)*(sin(pi*t)/(4*pi) + 1/4) + diffpotential*pi^2*sin(pi*y)*(sin(pi*t)/(4*pi) + 1/4) - 1'
  []

  [em_ICs]
    type = ParsedFunction
    vars = 'N_A'
    vals = 'N_A'
    value = 'log((3.0 + sin(pi/2*x)) / N_A)'
  []
  [ion_ICs]
    type = ParsedFunction
    vars = 'N_A'
    vals = 'N_A'
    value = 'log((3.0 + sin(pi/2*x)) / N_A)'
  []
  [mean_en_ICs]
    type = ParsedFunction
    vars = 'em_ICs'
    vals = 'em_ICs'
    value = 'log(1.0)'
  []

  [V_func_left_BC]
    type = ParsedFunction
    vars = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    vals = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    value = '-diffpotential*pi*cos(pi*x)*(sin(pi*t)/(4*pi) + 1/4)'
  []
  [V_func_down_BC]
    type = ParsedFunction
    vars = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    vals = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    value = '-diffpotential*pi*cos(pi*y)*(sin(pi*t)/(4*pi) + 1/4)'
  []
[]

[BCs]
  [em_BC]
    type = FunctionDirichletBC
    variable = em
    function = 'em_fun'
    boundary = '0 1 2 3'
    preset = true
  []

  [ion_BC]
    type = FunctionDirichletBC
    variable = ion
    function = 'ion_fun'
    boundary = '0 1 2 3'
    preset = true
  []

  [energy_BC]
    type = FunctionDirichletBC
    variable = mean_en
    function = 'mean_en_fun'
    boundary = '0 1 2 3'
    preset = true
  []

  [Ex_BC]
    type = FunctionDirichletBC
    variable = Ex
    function = 'Ex_fun'
    boundary = '0 1 2 3'
    preset = true
  []

  [Ey_BC]
    type = FunctionDirichletBC
    variable = Ey
    function = 'Ey_fun'
    boundary = '0 1 2 3'
    preset = true
  []

  [potential_left_BC]
    type = DielectricBCWithEffEfield
    variable = potential
    electrons = em
    electron_energy = mean_en
    ions = 'ion'
    Ex = Ex
    Ey = Ey
    dielectric_constant = 1.0
    thickness = 1.0
    emission_coeffs = 'users_gamma'
    potential_units = V
    position_units = 1.0
    boundary = 3
  []

  [potential_Dirichlet_BC]
    type = FunctionDirichletBC
    variable = potential
    function = 'potential_fun'
    boundary = '0 1 2'
    preset = true
  []
[]

[Materials]
  [field_solver]
    type = FieldSolverMaterial
    potential = potential
  []
  [Material_Coeff]
    type = GenericFunctionMaterial
    prop_names = 'e  N_A  massem NonAD_muion  NonAD_diffpotential  diffEx         diffEy'
    prop_values = 'ee N_A  massem muion        diffpotential        diffpotential  diffpotential '
  []
  [ADMaterial_Coeff_Set1]
    type = ADGenericFunctionMaterial
    prop_names = 'diffion  muion  diffem  muem  diffmean_en  diffpotential'
    prop_values = 'diffion  muion  diffem  muem  diffmean_en  diffpotential'
  []
  [Charge_Signs]
    type = GenericConstantMaterial
    prop_names = 'sgnem  sgnion  sgnmean_en'
    prop_values = '-1.0   1.0     -1.0'
  []
  [emission_coeffs]
    type = ADGenericConstantMaterial
    prop_names = 'users_gamma'
    prop_values = '1.0'
  []
[]

[Postprocessors]
  [em_l2Error]
    type = ElementL2Error
    variable = em
    function = em_fun
  []
  [ion_l2Error]
    type = ElementL2Error
    variable = ion
    function = ion_fun
  []
  [mean_en_l2Error]
    type = ElementL2Error
    variable = mean_en
    function = mean_en_fun
  []

  [Ex_l2Error]
    type = ElementL2Error
    variable = Ex
    function = Ex_fun
  []
  [Ey_l2Error]
    type = ElementL2Error
    variable = Ey
    function = Ey_fun
  []

  [potential_l2Error]
    type = ElementL2Error
    variable = potential
    function = potential_fun
  []

  [h]
    type = AverageElementSize
  []
[]

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

  [fdp]
    type = FDP
    full = true
  []
[]

[Executioner]
  type = Transient
  # end_time = 50
  start_time = 50
  end_time = 51

  # dt = 0.008
  dt = 0.02

  automatic_scaling = true
  compute_scaling_once = false
  petsc_options = '-snes_converged_reason -snes_linesearch_monitor'
  solve_type = NEWTON
  line_search = none
  petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount'
  petsc_options_value = 'lu NONZERO 1.e-10'

  scheme = bdf2

  nl_abs_tol = 1e-13
[]

[Outputs]
  [out]
    type = Exodus
    interval = 10
  []
[]

(test/tests/mms/bcs/2D_DielectricBCWithEffEfield.i)

[Mesh]
  [geo]
    type = FileMeshGenerator
    file = '2D_DielectricBCWithEffEfield_IC_out.e'
    use_for_exodus_restart = true
  []
[]

[Problem]
  type = FEProblem
[]

[Variables]
  [em]
    initial_from_file_var = em
  []
  [ion]
    initial_from_file_var = ion
  []
  [mean_en]
    initial_from_file_var = mean_en
  []

  [Ex]
    initial_from_file_var = Ex
  []
  [Ey]
    initial_from_file_var = Ey
  []

  [potential]
    initial_from_file_var = potential
  []
[]

# [ICs]
#   [em_IC]
#     type = FunctionIC
#     variable = em
#     function = 'em_ICs'
#   []
#   [ion_IC]
#     type = FunctionIC
#     variable = ion
#     function = 'ion_ICs'
#   []
#   [mean_en_IC]
#     type = FunctionIC
#     variable = mean_en
#     function = 'mean_en_ICs'
#   []
# []

[Kernels]
  #Electron Equations
  [em_time_derivative]
    type = ElectronTimeDerivative
    variable = em
  []
  [em_diffusion]
    type = CoeffDiffusion
    variable = em
    position_units = 1.0
  []
  [em_advection]
    type = EFieldAdvection
    variable = em
    position_units = 1.0
  []
  [em_source]
    type = BodyForce
    variable = em
    function = 'em_source'
  []

  #Ion Equations
  [ion_time_derivative]
    type = ElectronTimeDerivative
    variable = ion
  []
  [ion_diffusion]
    type = CoeffDiffusion
    variable = ion
    position_units = 1.0
  []
  [ion_advection]
    type = EffectiveEFieldAdvection
    variable = ion
    u = Ex
    v = Ey
    position_units = 1.0
  []
  [ion_source]
    type = BodyForce
    variable = ion
    function = 'ion_source'
  []

  #Eff. Efield
  [EffEfield_X_time_deriv]
    type = ADTimeDerivative
    variable = Ex
  []
  [EffEfield_X_diffusion]
    type = MatDiffusion
    D_name = diffEx
    variable = Ex
  []
  [EffEfield_X_source]
    type = BodyForce
    variable = Ex
    function = 'Ex_source'
  []
  [EffEfield_Y_time_deriv]
    type = ADTimeDerivative
    variable = Ey
  []
  [EffEfield_Y_diffusion]
    type = MatDiffusion
    D_name = diffEy
    variable = Ey
  []
  [EffEfield_Y_source]
    type = BodyForce
    variable = Ey
    function = 'Ey_source'
  []

  #Potential
  [Potential_time_deriv]
    type = ADTimeDerivative
    variable = potential
  []
  [Potential_diffusion]
    type = CoeffDiffusionLin
    variable = potential
    position_units = 1.0
  []
  [Potential_source]
    type = BodyForce
    variable = potential
    function = 'potential_source'
  []

  #Electron Energy Equations
  [mean_en_time_deriv]
    type = ElectronTimeDerivative
    variable = mean_en
  []
  [mean_en_diffusion]
    type = CoeffDiffusion
    variable = mean_en
    position_units = 1.0
  []
  [mean_en_source]
    type = BodyForce
    variable = mean_en
    function = 'energy_source'
  []
[]

[AuxVariables]
  [mean_en_sol]
  []

  [em_sol]
  []

  [ion_sol]
  []

  [Ex_sol]
  []
  [Ey_sol]
  []

  [potential_sol]
  []
[]

[AuxKernels]
  [mean_en_sol]
    type = FunctionAux
    variable = mean_en_sol
    function = mean_en_fun
  []

  [em_sol]
    type = FunctionAux
    variable = em_sol
    function = em_fun
  []

  [ion_sol]
    type = FunctionAux
    variable = ion_sol
    function = ion_fun
  []

  [Ex_sol]
    type = FunctionAux
    variable = Ex_sol
    function = Ex_fun
  []
  [Ey_sol]
    type = FunctionAux
    variable = Ey_sol
    function = Ey_fun
  []

  [potential_sol]
    type = FunctionAux
    variable = potential_sol
    function = potential_fun
  []
[]

[Functions]
  #Material Variables
  [massem]
    type = ConstantFunction
    value = 2.0
  []
  #Electron diffusion coeff.
  [diffem]
    type = ConstantFunction
    value = 0.05
  []
  [muem]
    type = ConstantFunction
    value = 0.01
  []
  #Electron energy mobility coeff.
  [diffmean_en]
    type = ConstantFunction
    value = 0.05
  []
  #Ion diffusion coeff.
  [diffion]
    type = ParsedFunction
    vars = diffem
    vals = diffem
    value = diffem
  []
  [muion]
    type = ParsedFunction
    vars = muem
    vals = muem
    value = muem
  []
  [N_A]
    type = ConstantFunction
    value = 1.0
  []
  [ee]
    type = ConstantFunction
    value = 1.0
  []
  [diffpotential]
    type = ConstantFunction
    value = 0.05
  []

  #Manufactured Solutions
  #The manufactured electron density solution
  [em_fun]
    type = ParsedFunction
    vars = 'N_A'
    vals = 'N_A'
    value = 'log((sin(pi*y) + 0.2*sin(2*pi*t)*cos(pi*y) + 2.0 + sin(pi*x)) / N_A)'
  []
  #The manufactured ion density solution
  [ion_fun]
    type = ParsedFunction
    vars = 'N_A'
    vals = 'N_A'
    value = 'log((sin(pi*y) + 0.2*sin(2*pi*t)*cos(pi*y) + 2.0 + sin(pi*x)) / N_A)'
  []
  #The manufactured electron energy solution
  [mean_en_fun]
    type = ParsedFunction
    vars = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    vals = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    value = 'log(((3*massem*pi*((diffpotential*pi*cos(pi*t) - cos(pi*t)*(sin(pi*x) + sin(pi*y)) + 4)/ee + 8*muion*pi*(sin(pi*t) + 1)*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2))^2)/(16*ee*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2))) / N_A)'
  []
  #The manufactured eff. Efield solution
  [Ex_fun]
    type = ParsedFunction
    vars = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    vals = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    value = '-pi*cos(pi*x)*(sin(pi*t) + 1)'
  []
  [Ey_fun]
    type = ParsedFunction
    vars = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    vals = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    value = '-pi*cos(pi*y)*(sin(pi*t) + 1)'
  []
  #The manufactured potential solution
  [potential_fun]
    type = ParsedFunction
    vars = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    vals = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    value = '0.25*(sin(pi*t)/pi + 1.0)*(sin(pi*y) + sin(pi*x)) - t'
  []

  #Source Terms in moles
  #The electron source term.
  [em_source]
    type = ParsedFunction
    vars = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    vals = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    value = '(diffem*pi^2*sin(pi*x) + (diffem*pi^2*(5*sin(pi*y) + cos(pi*y)*sin(2*pi*t)))/5 + (2*pi*cos(2*pi*t)*cos(pi*y))/5 - (muem*pi*(pi + sin(pi*t))*(10*sin(pi*x) + 10*sin(pi*y) + 10*sin(pi*x)*sin(pi*y) - 10*cos(pi*x)^2 - 10*cos(pi*y)^2 + cos(pi*y)*sin(2*pi*t)*sin(pi*x) + 2*cos(pi*y)*sin(2*pi*t)*sin(pi*y) + 10))/20) / N_A'
  []
  #The ion source term.
  [ion_source]
    type = ParsedFunction
    vars = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    vals = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    value = '(diffion*pi^2*sin(pi*x) + (diffion*pi^2*(5*sin(pi*y) + cos(pi*y)*sin(2*pi*t)))/5 + (2*pi*cos(2*pi*t)*cos(pi*y))/5 + (muion*pi^2*(sin(pi*t) + 1)*(10*sin(pi*x) + 10*sin(pi*y) + 10*sin(pi*x)*sin(pi*y) - 10*cos(pi*x)^2 - 10*cos(pi*y)^2 + cos(pi*y)*sin(2*pi*t)*sin(pi*x) + 2*cos(pi*y)*sin(2*pi*t)*sin(pi*y) + 10))/5) / N_A'
  []
  [energy_source]
    type = ParsedFunction
    vars = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    vals = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    value = '((3*massem*pi*((diffpotential*pi*cos(pi*t) - cos(pi*t)*(sin(pi*x) + sin(pi*y)) + 4)/ee + 8*muion*pi*(sin(pi*t) + 1)*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2))*(8*muion*pi^2*cos(pi*t)*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2) - (diffpotential*pi^2*sin(pi*t) - pi*sin(pi*t)*(sin(pi*x) + sin(pi*y)))/ee + (16*muion*pi^2*cos(2*pi*t)*cos(pi*y)*(sin(pi*t) + 1))/5))/(8*ee*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2)) - diffmean_en*((3*massem*pi*((pi*cos(pi*t)*cos(pi*y))/ee - 8*muion*pi*(pi*cos(pi*y) - (pi*sin(2*pi*t)*sin(pi*y))/5)*(sin(pi*t) + 1))^2)/(8*ee*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2)) + (3*massem*pi*((diffpotential*pi*cos(pi*t) - cos(pi*t)*(sin(pi*x) + sin(pi*y)) + 4)/ee + 8*muion*pi*(sin(pi*t) + 1)*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2))^2*(pi^2*sin(pi*y) + (pi^2*cos(pi*y)*sin(2*pi*t))/5))/(16*ee*(sin(x*pi) + sin(y*pi) + (cos(y*pi)*sin(2*t*pi))/5 + 2)^2) - (3*massem*pi*((diffpotential*pi*cos(pi*t) - cos(pi*t)*(sin(pi*x) + sin(pi*y)) + 4)/ee + 8*muion*pi*(sin(pi*t) + 1)*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2))*(8*muion*pi*(pi^2*sin(pi*y) + (pi^2*cos(pi*y)*sin(2*pi*t))/5)*(sin(pi*t) + 1) - (pi^2*cos(pi*t)*sin(pi*y))/ee))/(8*ee*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2)) + (3*massem*pi*((diffpotential*pi*cos(pi*t) - cos(pi*t)*(sin(pi*x) + sin(pi*y)) + 4)/ee + 8*muion*pi*(sin(pi*t) + 1)*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2))^2*(pi*cos(pi*y) - (pi*sin(2*pi*t)*sin(pi*y))/5)^2)/(8*ee*(sin(x*pi) + sin(y*pi) + (cos(y*pi)*sin(2*t*pi))/5 + 2)^3) + (3*massem*pi*((diffpotential*pi*cos(pi*t) - cos(pi*t)*(sin(pi*x) + sin(pi*y)) + 4)/ee + 8*muion*pi*(sin(pi*t) + 1)*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2))*((pi*cos(pi*t)*cos(pi*y))/ee - 8*muion*pi*(pi*cos(pi*y) - (pi*sin(2*pi*t)*sin(pi*y))/5)*(sin(pi*t) + 1))*(pi*cos(pi*y) - (pi*sin(2*pi*t)*sin(pi*y))/5))/(4*ee*(sin(x*pi) + sin(y*pi) + (cos(y*pi)*sin(2*t*pi))/5 + 2)^2)) - diffmean_en*((3*massem*pi*((pi*cos(pi*t)*cos(pi*x))/ee - 8*muion*pi^2*cos(pi*x)*(sin(pi*t) + 1))^2)/(8*ee*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2)) + (3*massem*pi^3*cos(pi*x)^2*((diffpotential*pi*cos(pi*t) - cos(pi*t)*(sin(pi*x) + sin(pi*y)) + 4)/ee + 8*muion*pi*(sin(pi*t) + 1)*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2))^2)/(8*ee*(sin(x*pi) + sin(y*pi) + (cos(y*pi)*sin(2*t*pi))/5 + 2)^3) + (3*massem*pi^3*sin(pi*x)*((diffpotential*pi*cos(pi*t) - cos(pi*t)*(sin(pi*x) + sin(pi*y)) + 4)/ee + 8*muion*pi*(sin(pi*t) + 1)*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2))^2)/(16*ee*(sin(x*pi) + sin(y*pi) + (cos(y*pi)*sin(2*t*pi))/5 + 2)^2) - (3*massem*pi*((diffpotential*pi*cos(pi*t) - cos(pi*t)*(sin(pi*x) + sin(pi*y)) + 4)/ee + 8*muion*pi*(sin(pi*t) + 1)*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2))*(8*muion*pi^3*sin(pi*x)*(sin(pi*t) + 1) - (pi^2*cos(pi*t)*sin(pi*x))/ee))/(8*ee*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2)) + (3*massem*pi^2*cos(pi*x)*((diffpotential*pi*cos(pi*t) - cos(pi*t)*(sin(pi*x) + sin(pi*y)) + 4)/ee + 8*muion*pi*(sin(pi*t) + 1)*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2))*((pi*cos(pi*t)*cos(pi*x))/ee - 8*muion*pi^2*cos(pi*x)*(sin(pi*t) + 1)))/(4*ee*(sin(x*pi) + sin(y*pi) + (cos(y*pi)*sin(2*t*pi))/5 + 2)^2)) - (3*massem*pi^2*cos(2*pi*t)*cos(pi*y)*((diffpotential*pi*cos(pi*t) - cos(pi*t)*(sin(pi*x) + sin(pi*y)) + 4)/ee + 8*muion*pi*(sin(pi*t) + 1)*(sin(pi*x) + sin(pi*y) + (cos(pi*y)*sin(2*pi*t))/5 + 2))^2)/(40*ee*(sin(x*pi) + sin(y*pi) + (cos(y*pi)*sin(2*t*pi))/5 + 2)^2)) / N_A'
  []

  #The Ex source term.
  [Ex_source]
    type = ParsedFunction
    vars = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    vals = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    value = '-pi^2*cos(pi*t)*cos(pi*x) - diffpotential*pi^3*cos(pi*x)*(sin(pi*t) + 1)'
  []
  [Ey_source]
    type = ParsedFunction
    vars = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    vals = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    value = '-pi^2*cos(pi*t)*cos(pi*y) - diffpotential*pi^3*cos(pi*y)*(sin(pi*t) + 1)'
  []

  [potential_source]
    type = ParsedFunction
    vars = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    vals = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    value = '(cos(pi*t)*(sin(pi*x) + sin(pi*y)))/4 + diffpotential*pi^2*sin(pi*x)*(sin(pi*t)/(4*pi) + 1/4) + diffpotential*pi^2*sin(pi*y)*(sin(pi*t)/(4*pi) + 1/4) - 1'
  []

  [em_ICs]
    type = ParsedFunction
    vars = 'N_A'
    vals = 'N_A'
    value = 'log((3.0 + sin(pi/2*x)) / N_A)'
  []
  [ion_ICs]
    type = ParsedFunction
    vars = 'N_A'
    vals = 'N_A'
    value = 'log((3.0 + sin(pi/2*x)) / N_A)'
  []
  [mean_en_ICs]
    type = ParsedFunction
    vars = 'em_ICs'
    vals = 'em_ICs'
    value = 'log(1.0)'
  []

  [V_func_left_BC]
    type = ParsedFunction
    vars = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    vals = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    value = '-diffpotential*pi*cos(pi*x)*(sin(pi*t)/(4*pi) + 1/4)'
  []
  [V_func_down_BC]
    type = ParsedFunction
    vars = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    vals = 'ee N_A diffpotential diffem muem massem diffmean_en diffion muion'
    value = '-diffpotential*pi*cos(pi*y)*(sin(pi*t)/(4*pi) + 1/4)'
  []
[]

[BCs]
  [em_BC]
    type = FunctionDirichletBC
    variable = em
    function = 'em_fun'
    boundary = '0 1 2 3'
    preset = true
  []

  [ion_BC]
    type = FunctionDirichletBC
    variable = ion
    function = 'ion_fun'
    boundary = '0 1 2 3'
    preset = true
  []

  [energy_BC]
    type = FunctionDirichletBC
    variable = mean_en
    function = 'mean_en_fun'
    boundary = '0 1 2 3'
    preset = true
  []

  [Ex_BC]
    type = FunctionDirichletBC
    variable = Ex
    function = 'Ex_fun'
    boundary = '0 1 2 3'
    preset = true
  []

  [Ey_BC]
    type = FunctionDirichletBC
    variable = Ey
    function = 'Ey_fun'
    boundary = '0 1 2 3'
    preset = true
  []

  [potential_left_BC]
    type = DielectricBCWithEffEfield
    variable = potential
    electrons = em
    electron_energy = mean_en
    ions = 'ion'
    Ex = Ex
    Ey = Ey
    dielectric_constant = 1.0
    thickness = 1.0
    emission_coeffs = 'users_gamma'
    potential_units = V
    position_units = 1.0
    boundary = 3
  []

  [potential_Dirichlet_BC]
    type = FunctionDirichletBC
    variable = potential
    function = 'potential_fun'
    boundary = '0 1 2'
    preset = true
  []
[]

[Materials]
  [field_solver]
    type = FieldSolverMaterial
    potential = potential
  []
  [Material_Coeff]
    type = GenericFunctionMaterial
    prop_names = 'e  N_A  massem NonAD_muion  NonAD_diffpotential  diffEx         diffEy'
    prop_values = 'ee N_A  massem muion        diffpotential        diffpotential  diffpotential '
  []
  [ADMaterial_Coeff_Set1]
    type = ADGenericFunctionMaterial
    prop_names = 'diffion  muion  diffem  muem  diffmean_en  diffpotential'
    prop_values = 'diffion  muion  diffem  muem  diffmean_en  diffpotential'
  []
  [Charge_Signs]
    type = GenericConstantMaterial
    prop_names = 'sgnem  sgnion  sgnmean_en'
    prop_values = '-1.0   1.0     -1.0'
  []
  [emission_coeffs]
    type = ADGenericConstantMaterial
    prop_names = 'users_gamma'
    prop_values = '1.0'
  []
[]

[Postprocessors]
  [em_l2Error]
    type = ElementL2Error
    variable = em
    function = em_fun
  []
  [ion_l2Error]
    type = ElementL2Error
    variable = ion
    function = ion_fun
  []
  [mean_en_l2Error]
    type = ElementL2Error
    variable = mean_en
    function = mean_en_fun
  []

  [Ex_l2Error]
    type = ElementL2Error
    variable = Ex
    function = Ex_fun
  []
  [Ey_l2Error]
    type = ElementL2Error
    variable = Ey
    function = Ey_fun
  []

  [potential_l2Error]
    type = ElementL2Error
    variable = potential
    function = potential_fun
  []

  [h]
    type = AverageElementSize
  []
[]

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

  [fdp]
    type = FDP
    full = true
  []
[]

[Executioner]
  type = Transient
  # end_time = 50
  start_time = 50
  end_time = 51

  # dt = 0.008
  dt = 0.02

  automatic_scaling = true
  compute_scaling_once = false
  petsc_options = '-snes_converged_reason -snes_linesearch_monitor'
  solve_type = NEWTON
  line_search = none
  petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount'
  petsc_options_value = 'lu NONZERO 1.e-10'

  scheme = bdf2

  nl_abs_tol = 1e-13
[]

[Outputs]
  [out]
    type = Exodus
    interval = 10
  []
[]

Overview
Example Input File Syntax
Input Parameters
Input Files
References