ThermalContactCondition

Interface condition that describes the thermal contact resistance across a boundary formed between two dissimilar materials (resulting in a temperature discontinuity) under the influence of an electrostatic potential, as described in Cincotti, et al (DOI: 10.1002/aic.11102). Thermal conductivity on each side of the boundary is defined via the material properties system.

Description

This interface kernel models the conductivity of heat flux across a specified boundary between two dissimilar materials, as described by (Cincotti et al., 2007). It accounts for the influence of both temperature and electrostatic potential differences across the interface, with appropriate thermal and electrical contact conductances being provided by the user as a constant scalar number or via a combination of material properties and constants for calculation. The condition being applied is:

$var element = document.getElementById("moose-equation-24dbf684-33b8-4d94-83bb-171e805882b3");katex.render(" -k_1 \\frac{\\partial T}{\\partial \\mathbf{r}} \\bigg\\rvert_1 \\cdot \\mathbf{\\hat{n}} = C_T (T_1 - T_2) - f C_E (\\phi_1 - \\phi_2)^2", element, {displayMode:true,throwOnError:false});$

$var element = document.getElementById("moose-equation-7fd3dc61-1711-4447-b92e-bc9e11c15b64");katex.render(" -k_2 \\frac{\\partial T}{\\partial \\mathbf{r}} \\bigg\\rvert_2 \\cdot \\mathbf{\\hat{n}} = C_T (T_1 - T_2) + (1 - f) C_E (\\phi_1 - \\phi_2)^2", element, {displayMode:true,throwOnError:false});$

where

$var element = document.getElementById("moose-equation-fadd3cf7-6744-4508-9cfb-16e3e980aea7");katex.render("k_i", element, {displayMode:false,throwOnError:false});$ is the thermal conductivity of each material along the interface,
$var element = document.getElementById("moose-equation-a8c568b2-70b7-497c-b18e-0462493b595f");katex.render("C_T", element, {displayMode:false,throwOnError:false});$ is the thermal contact conductance,
$var element = document.getElementById("moose-equation-e5c6e0fe-b952-409e-b7ce-0dc7a2a81ca9");katex.render("C_E", element, {displayMode:false,throwOnError:false});$ is the electrical contact conductance,
$var element = document.getElementById("moose-equation-dc35a2f2-bdb2-47bd-aa0c-7987946bb2d4");katex.render("T_i", element, {displayMode:false,throwOnError:false});$ is the temperature of the material at the interface,
$var element = document.getElementById("moose-equation-108f83f9-8f5c-4580-a0d0-cfc339b8a0f1");katex.render("\\phi_i", element, {displayMode:false,throwOnError:false});$ is the electrostatic potential of the material at the interface, and
$var element = document.getElementById("moose-equation-18ac03c3-e543-4f20-a22e-d85733ff1e5f");katex.render("f", element, {displayMode:false,throwOnError:false});$ is the splitting factor of the Joule heat source arising from electrical contact resistance (default is 0.5). This parameter sets the fraction of heat flux associated with Joule heating that enters the primary side of the boundary, with the rest entering the secondary side.

Thermal Contact Conductance

The temperature- and mechanical-pressure-dependent thermal contact conductance, given by (Madhusudana, 1996), is calculated using:

$var element = document.getElementById("moose-equation-1c9e9e29-ad29-490f-a071-0b6404873f18");katex.render(" C_T(T, P) = \\alpha_T k_{Harm} \\bigg( \\frac{P}{H_{Harm}} \\bigg)^{\\beta_T}", element, {displayMode:true,throwOnError:false});$

where

$var element = document.getElementById("moose-equation-befcd875-f17d-4a08-965a-00828d1c8ef0");katex.render("\\alpha_T", element, {displayMode:false,throwOnError:false});$ is an experimentally-derived proportional fit parameter (set to be $var element = document.getElementById("moose-equation-88adfff5-5ab5-4f18-aedd-1645c8bfb112");katex.render("22,810", element, {displayMode:false,throwOnError:false});$ , from (Cincotti et al., 2007)),
$var element = document.getElementById("moose-equation-49be45b5-3447-484f-8e4d-6a7910f04355");katex.render("k_{Harm}", element, {displayMode:false,throwOnError:false});$ is the harmonic mean of the temperature-dependent thermal conductivities on either side of the boundary,
$var element = document.getElementById("moose-equation-ea2c90c3-111d-4075-b97a-804737171da5");katex.render("P", element, {displayMode:false,throwOnError:false});$ ( $var element = document.getElementById("moose-equation-f1ed2a24-6e6c-449f-8d89-96d4b196a5b9");katex.render("=F/S", element, {displayMode:false,throwOnError:false});$ ) is the uniform mechanical pressure applied at the contact surface area (S) between the two materials,
$var element = document.getElementById("moose-equation-da7e4525-eaf0-4a1f-888c-a145e090af11");katex.render("H_{Harm}", element, {displayMode:false,throwOnError:false});$ is the harmonic mean of the hardness values of each material, and
$var element = document.getElementById("moose-equation-2e857c3d-6cf9-438f-9411-56e2c486ff46");katex.render("\\beta_T", element, {displayMode:false,throwOnError:false});$ is an experimentally-derived power fit parameter (set to be $var element = document.getElementById("moose-equation-6211f1c7-3b92-42a6-a0b0-d4c4340f30dd");katex.render("1.08", element, {displayMode:false,throwOnError:false});$ , from (Cincotti et al., 2007)).

Electrical Contact Conductance

The temperature- and mechanical-pressure-dependent electrical contact conductance, given by (Babu et al., 2001), is calculated using:

$var element = document.getElementById("moose-equation-6bc8e212-0544-4ffc-ac56-e801fb60a66f");katex.render(" C_E(T, P) = \\alpha_E \\sigma_{el,Harm} \\bigg( \\frac{P}{H_{Harm}} \\bigg)^{\\beta_E}", element, {displayMode:true,throwOnError:false});$

where

$var element = document.getElementById("moose-equation-00a2c66a-2b00-4895-90e6-72bc13797934");katex.render("\\alpha_E", element, {displayMode:false,throwOnError:false});$ is an experimentally-derived proportional fit parameter (set to be $var element = document.getElementById("moose-equation-86241307-d4a5-43e9-970d-f3fc03eab6fe");katex.render("64", element, {displayMode:false,throwOnError:false});$ , from (Cincotti et al., 2007)),
$var element = document.getElementById("moose-equation-a69503af-dc57-49f3-9fc7-992566bcf4b5");katex.render("\\sigma_{el,Harm}", element, {displayMode:false,throwOnError:false});$ is the harmonic mean of the temperature-dependent electrical conductivities on either side of the boundary,
$var element = document.getElementById("moose-equation-0d61e97d-7a84-40a9-98d7-d8a0828082bc");katex.render("P", element, {displayMode:false,throwOnError:false});$ ( $var element = document.getElementById("moose-equation-99c70276-eaef-42ee-98e4-c1ecd26db6a0");katex.render("=F/S", element, {displayMode:false,throwOnError:false});$ ) is the uniform mechanical pressure applied at the contact surface area (S) between the two materials,
$var element = document.getElementById("moose-equation-2e23fba2-06bc-4a5d-8d38-68b33be9de8c");katex.render("H_{Harm}", element, {displayMode:false,throwOnError:false});$ is the harmonic mean of the hardness values of each material, and
$var element = document.getElementById("moose-equation-bce48950-ad8f-4473-a735-44b2793625b6");katex.render("\\beta_E", element, {displayMode:false,throwOnError:false});$ is an experimentally-derived power fit parameter (set to be $var element = document.getElementById("moose-equation-d1836efa-b5b1-47bd-96ac-8e7effb31a04");katex.render("0.35", element, {displayMode:false,throwOnError:false});$ , from (Cincotti et al., 2007)).

Geometric Mean

For reference, the harmonic mean calculation for two values, $var element = document.getElementById("moose-equation-cb7a9095-e20b-43ee-a8cc-fbf9a0525160");katex.render("V_a", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-5fe5e897-2537-45a2-b4be-c3127fcddd18");katex.render("V_b", element, {displayMode:false,throwOnError:false});$ , is given by:

$var element = document.getElementById("moose-equation-4076c2cb-1555-4951-beff-a4179e4e448e");katex.render(" V_{Harm} = \\frac{2 V_a V_b}{V_a + V_b}", element, {displayMode:true,throwOnError:false});$

Example Input File Syntax

Calculated Conductance

[InterfaceKernels<<<{"href": "../../syntax/InterfaceKernels/index.html"}>>>]
  [thermal_contact_conductance_calculated]
    type = ThermalContactCondition<<<{"description": "Interface condition that describes the thermal contact resistance across a boundary formed between two dissimilar materials (resulting in a temperature discontinuity) under the influence of an electrostatic potential, as described in Cincotti, et al (DOI: 10.1002/aic.11102). Thermal conductivity on each side of the boundary is defined via the material properties system.", "href": "ThermalContactCondition.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature_stainless_steel
    neighbor_var<<<{"description": "The variable on the other side of the interface."}>>> = temperature_graphite
    primary_potential<<<{"description": "Electrostatic potential on the primary block."}>>> = potential_stainless_steel
    secondary_potential<<<{"description": "Electrostatic potential on the secondary block."}>>> = potential_graphite
    primary_thermal_conductivity<<<{"description": "Thermal conductivity on the primary block."}>>> = thermal_conductivity
    secondary_thermal_conductivity<<<{"description": "Thermal conductivity on the secondary block."}>>> = thermal_conductivity
    primary_electrical_conductivity<<<{"description": "Electrical conductivity on the primary block."}>>> = electrical_conductivity
    secondary_electrical_conductivity<<<{"description": "Electrical conductivity on the secondary block."}>>> = electrical_conductivity
    mean_hardness<<<{"description": "Geometric mean of the hardness of each contacting material."}>>> = graphite_stainless_mean_hardness
    mechanical_pressure<<<{"description": "Mechanical pressure uniformly applied at the contact surface area (Pressure = Force / Surface Area)."}>>> = 8.52842e10 # resulting in electrical contact conductance = ~1.4715e5, thermal contact conductance = ~3.44689e7
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = ssg_interface
  []
[]

Supplied Conductance

[InterfaceKernels<<<{"href": "../../syntax/InterfaceKernels/index.html"}>>>]
  [thermal_contact_conductance]
    type = ThermalContactCondition<<<{"description": "Interface condition that describes the thermal contact resistance across a boundary formed between two dissimilar materials (resulting in a temperature discontinuity) under the influence of an electrostatic potential, as described in Cincotti, et al (DOI: 10.1002/aic.11102). Thermal conductivity on each side of the boundary is defined via the material properties system.", "href": "ThermalContactCondition.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature_stainless_steel
    neighbor_var<<<{"description": "The variable on the other side of the interface."}>>> = temperature_graphite
    primary_potential<<<{"description": "Electrostatic potential on the primary block."}>>> = potential_stainless_steel
    secondary_potential<<<{"description": "Electrostatic potential on the secondary block."}>>> = potential_graphite
    primary_thermal_conductivity<<<{"description": "Thermal conductivity on the primary block."}>>> = thermal_conductivity
    secondary_thermal_conductivity<<<{"description": "Thermal conductivity on the secondary block."}>>> = thermal_conductivity
    primary_electrical_conductivity<<<{"description": "Electrical conductivity on the primary block."}>>> = electrical_conductivity
    secondary_electrical_conductivity<<<{"description": "Electrical conductivity on the secondary block."}>>> = electrical_conductivity
    user_electrical_contact_conductance<<<{"description": "User-provided electrical contact conductance coefficient."}>>> = 2.5e5 # as described in Cincotti et al (DOI: 10.1002/aic.11102)
    user_thermal_contact_conductance<<<{"description": "User-provided thermal contact conductance coefficient."}>>> = 7 # also from Cincotti et al
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = ssg_interface
  []
[]

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.
primary_potentialElectrostatic potential on the primary block.
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Electrostatic potential on the primary block.
secondary_potentialElectrostatic potential on the secondary block.
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Electrostatic potential on the secondary block.
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)
mean_hardnessmean_hardnessGeometric mean of the hardness of each contacting material.
Default:mean_hardness
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:Geometric mean of the hardness of each contacting material.
mechanical_pressure0Mechanical pressure uniformly applied at the contact surface area (Pressure = Force / Surface Area).
Default:0
C++ Type:FunctionName
Unit:(no unit assumed)
Controllable:No
Description:Mechanical pressure uniformly applied at the contact surface area (Pressure = Force / Surface Area).
primary_electrical_conductivityelectrical_conductivityElectrical conductivity on the primary block.
Default:electrical_conductivity
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:Electrical conductivity on the primary block.
primary_thermal_conductivitythermal_conductivityThermal conductivity on the primary block.
Default:thermal_conductivity
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:Thermal conductivity on the primary block.
secondary_electrical_conductivityelectrical_conductivityElectrical conductivity on the secondary block.
Default:electrical_conductivity
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:Electrical conductivity on the secondary block.
secondary_thermal_conductivitythermal_conductivityThermal conductivity on the secondary block.
Default:thermal_conductivity
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:Thermal conductivity on the secondary block.
splitting_factor0.5Splitting factor of the Joule heating source.
Default:0.5
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Splitting factor of the Joule heating source.
user_electrical_contact_conductance0User-provided electrical contact conductance coefficient.
Default:0
C++ Type:FunctionName
Unit:(no unit assumed)
Controllable:No
Description:User-provided electrical contact conductance coefficient.
user_thermal_contact_conductance0User-provided thermal contact conductance coefficient.
Default:0
C++ Type:FunctionName
Unit:(no unit assumed)
Controllable:No
Description:User-provided thermal contact conductance coefficient.

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/interfacekernels/thermal_conductance/thermal_interface.i)
(test/tests/example_testing/sample_three.i)
(examples/sps/multiapp/complete_geometry/electrothermal/engineering_scale_electrothermal.i)
(test/tests/example_testing/sample_four.i)
(test/tests/example_testing/sample_one.i)
(examples/sps/multiapp/electrothermal_with_phase_field/twoway_lots_of_particles_prototype/engineering_scale_electrothermal_twoway_lots_prototype.i)
(examples/sps/multiapp/complete_geometry/electrothermal/engineering_scale_electrothermal_goofycontact.i)
(test/tests/interfacekernels/thermal_conductance/thermal_interface_analytic_solution_two_block.i)
(examples/sps/multiapp/electrothermal_with_phase_field/twoway_initial_prototype/engineering_scale_electrothermal_twoway_prototype.i)
(examples/sps/multiapp/electrothermal_with_phase_field/oneway_controls/engineering_scale_electrothermal_oneway_controls.i)

References

S. S. Babu, M. L. Santella, Z. Feng, B. W. Riemer, and J. W. Cohron. Empirical model of effects of pressure and temperature on electrical contact resistance of metals. Sci Technol Weld Joining, 6(3):126–132, 2001.

@article{babu2001contactresistance,
    author = "Babu, S. S. and Santella, M. L. and Feng, Z. and Riemer, B. W. and Cohron, J. W.",
    title = "Empirical model of effects of pressure and temperature on electrical contact resistance of metals",
    journal = "Sci Technol Weld Joining",
    volume = "6",
    number = "3",
    pages = "126-132",
    year = "2001"
}

A. Cincotti, A. M. Locci, R. Orrù, and G. Cao. Modeling of SPS apparatus: temperature, current and strain distribution with no powders. AIChE Journal, 53(3):703–719, 2007. doi:10.1002/aic.11102.

@article{cincotti2007sps,
    author = "Cincotti, A. and Locci, A. M. and Orrù, R. and Cao, G.",
    title = "Modeling of {SPS} apparatus: Temperature, current and strain distribution with no powders",
    journal = "AIChE Journal",
    volume = "53",
    number = "3",
    pages = "703-719",
    doi = "10.1002/aic.11102",
    year = "2007"
}

Chakravarti V. Madhusudana. Thermal Contact Conductance. Springer-Verlag, 1996.

(test/tests/interfacekernels/thermal_conductance/thermal_interface.i)

# Input file block to generate this mesh can be found in thermal_interface_mesh.i
[Mesh]
  file = thermal_interface_regular_mesh.e
  coord_type = RZ
[]

[Variables]
  [temperature_graphite]
    initial_condition = 300.0
    block = graphite
  []
  [temperature_stainless_steel]
    initial_condition = 300.0
    block = stainless_steel
  []
  [potential_graphite]
    block = graphite
  []
  [potential_stainless_steel]
    block = stainless_steel
  []
[]

[AuxVariables]
  [T_infinity]
    initial_condition = 300.0
  []
[]

[Kernels]
  [HeatDiff_graphite]
    type = ADHeatConduction
    variable = temperature_graphite
    block = graphite
  []
  [HeatTdot_graphite]
    type = ADHeatConductionTimeDerivative
    variable = temperature_graphite
    block = graphite
  []
  [HeatSource_graphite]
    type = ADJouleHeatingSource
    variable = temperature_graphite
    elec = potential_graphite
    electrical_conductivity = electrical_conductivity
    block = graphite
  []

  [HeatDiff_stainless_steel]
    type = ADHeatConduction
    variable = temperature_stainless_steel
    block = stainless_steel
  []
  [HeatTdot_stainless_steel]
    type = ADHeatConductionTimeDerivative
    variable = temperature_stainless_steel
    block = stainless_steel
  []
  [HeatSource_stainless_steel]
    type = ADJouleHeatingSource
    variable = temperature_stainless_steel
    elec = potential_stainless_steel
    electrical_conductivity = electrical_conductivity
    block = stainless_steel
  []

  [electric_graphite]
    type = ADMatDiffusion
    variable = potential_graphite
    diffusivity = electrical_conductivity
    block = graphite
  []
  [electric_stainless_steel]
    type = ADMatDiffusion
    variable = potential_stainless_steel
    diffusivity = electrical_conductivity
    block = stainless_steel
  []
[]

[BCs]
  [external_surface_stainless]
    type = ADCoupledSimpleRadiativeHeatFluxBC
    boundary = right_stainless_steel
    variable = temperature_stainless_steel
    T_infinity = T_infinity
    emissivity = 0.85
  []
  [external_surface_graphite]
    type = ADCoupledSimpleRadiativeHeatFluxBC
    boundary = right_graphite
    variable = temperature_graphite
    T_infinity = T_infinity
    emissivity = 0.85
  []
  [elec_top]
    type = ADDirichletBC
    variable = potential_stainless_steel
    boundary = top_die
    value = 0.68 # Better reflects Cincotti et al (DOI: 10.1002/aic.11102) Figure 19
  []
  [elec_bottom]
    type = ADDirichletBC
    variable = potential_stainless_steel
    boundary = bottom_die
    value = 0
  []
[]

[InterfaceKernels]
  active = 'thermal_contact_conductance electric_contact_conductance'

  [thermal_contact_conductance]
    type = ThermalContactCondition
    variable = temperature_stainless_steel
    neighbor_var = temperature_graphite
    primary_potential = potential_stainless_steel
    secondary_potential = potential_graphite
    primary_thermal_conductivity = thermal_conductivity
    secondary_thermal_conductivity = thermal_conductivity
    primary_electrical_conductivity = electrical_conductivity
    secondary_electrical_conductivity = electrical_conductivity
    user_electrical_contact_conductance = 2.5e5 # as described in Cincotti et al (DOI: 10.1002/aic.11102)
    user_thermal_contact_conductance = 7 # also from Cincotti et al
    boundary = ssg_interface
  []
  [electric_contact_conductance]
    type = ElectrostaticContactCondition
    variable = potential_stainless_steel
    neighbor_var = potential_graphite
    primary_conductivity = electrical_conductivity
    secondary_conductivity = electrical_conductivity
    boundary = ssg_interface
    user_electrical_contact_conductance = 2.5e5 # as described in Cincotti et al (DOI: 10.1002/aic.11102)
  []

  [thermal_contact_conductance_calculated]
    type = ThermalContactCondition
    variable = temperature_stainless_steel
    neighbor_var = temperature_graphite
    primary_potential = potential_stainless_steel
    secondary_potential = potential_graphite
    primary_thermal_conductivity = thermal_conductivity
    secondary_thermal_conductivity = thermal_conductivity
    primary_electrical_conductivity = electrical_conductivity
    secondary_electrical_conductivity = electrical_conductivity
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = 8.52842e10 # resulting in electrical contact conductance = ~1.4715e5, thermal contact conductance = ~3.44689e7
    boundary = ssg_interface
  []
[]

[Materials]
  active = 'heat_conductor_graphite rho_graphite sigma_graphite heat_conductor_stainless_steel rho_stainless_steel sigma_stainless_steel'

  #graphite
  [heat_conductor_graphite]
    type = ADHeatConductionMaterial
    thermal_conductivity = 630
    specific_heat = 60
    block = graphite
  []
  [rho_graphite]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = 1.75e3
    block = graphite
  []
  [sigma_graphite]
    type = ADElectricalConductivity
    temperature = temperature_graphite
    reference_temperature = 293.0
    reference_resistivity = 3.0e-3
    temperature_coefficient = 0 # makes conductivity constant
    block = graphite
  []

  #stainless_steel
  [heat_conductor_stainless_steel]
    type = ADHeatConductionMaterial
    thermal_conductivity = 17
    specific_heat = 502
    block = stainless_steel
  []
  [rho_stainless_steel]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = 8e3
    block = stainless_steel
  []
  [sigma_stainless_steel]
    type = ADElectricalConductivity
    temperature = temperature_stainless_steel
    reference_temperature = 293.0
    reference_resistivity = 7e-7
    temperature_coefficient = 0 # makes conductivity constant
    block = stainless_steel
  []

  # harmonic mean of graphite and stainless steel hardness
  [mean_hardness]
    type = GraphiteStainlessMeanHardness
  []
[]

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

[Executioner]
  type = Transient
  scheme = bdf2
  solve_type = NEWTON
  dt = 0.1
  end_time = 100
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
[]

[Outputs]
  exodus = true
  perf_graph = true
[]

(test/tests/interfacekernels/thermal_conductance/thermal_interface.i)

# Input file block to generate this mesh can be found in thermal_interface_mesh.i
[Mesh]
  file = thermal_interface_regular_mesh.e
  coord_type = RZ
[]

[Variables]
  [temperature_graphite]
    initial_condition = 300.0
    block = graphite
  []
  [temperature_stainless_steel]
    initial_condition = 300.0
    block = stainless_steel
  []
  [potential_graphite]
    block = graphite
  []
  [potential_stainless_steel]
    block = stainless_steel
  []
[]

[AuxVariables]
  [T_infinity]
    initial_condition = 300.0
  []
[]

[Kernels]
  [HeatDiff_graphite]
    type = ADHeatConduction
    variable = temperature_graphite
    block = graphite
  []
  [HeatTdot_graphite]
    type = ADHeatConductionTimeDerivative
    variable = temperature_graphite
    block = graphite
  []
  [HeatSource_graphite]
    type = ADJouleHeatingSource
    variable = temperature_graphite
    elec = potential_graphite
    electrical_conductivity = electrical_conductivity
    block = graphite
  []

  [HeatDiff_stainless_steel]
    type = ADHeatConduction
    variable = temperature_stainless_steel
    block = stainless_steel
  []
  [HeatTdot_stainless_steel]
    type = ADHeatConductionTimeDerivative
    variable = temperature_stainless_steel
    block = stainless_steel
  []
  [HeatSource_stainless_steel]
    type = ADJouleHeatingSource
    variable = temperature_stainless_steel
    elec = potential_stainless_steel
    electrical_conductivity = electrical_conductivity
    block = stainless_steel
  []

  [electric_graphite]
    type = ADMatDiffusion
    variable = potential_graphite
    diffusivity = electrical_conductivity
    block = graphite
  []
  [electric_stainless_steel]
    type = ADMatDiffusion
    variable = potential_stainless_steel
    diffusivity = electrical_conductivity
    block = stainless_steel
  []
[]

[BCs]
  [external_surface_stainless]
    type = ADCoupledSimpleRadiativeHeatFluxBC
    boundary = right_stainless_steel
    variable = temperature_stainless_steel
    T_infinity = T_infinity
    emissivity = 0.85
  []
  [external_surface_graphite]
    type = ADCoupledSimpleRadiativeHeatFluxBC
    boundary = right_graphite
    variable = temperature_graphite
    T_infinity = T_infinity
    emissivity = 0.85
  []
  [elec_top]
    type = ADDirichletBC
    variable = potential_stainless_steel
    boundary = top_die
    value = 0.68 # Better reflects Cincotti et al (DOI: 10.1002/aic.11102) Figure 19
  []
  [elec_bottom]
    type = ADDirichletBC
    variable = potential_stainless_steel
    boundary = bottom_die
    value = 0
  []
[]

[InterfaceKernels]
  active = 'thermal_contact_conductance electric_contact_conductance'

  [thermal_contact_conductance]
    type = ThermalContactCondition
    variable = temperature_stainless_steel
    neighbor_var = temperature_graphite
    primary_potential = potential_stainless_steel
    secondary_potential = potential_graphite
    primary_thermal_conductivity = thermal_conductivity
    secondary_thermal_conductivity = thermal_conductivity
    primary_electrical_conductivity = electrical_conductivity
    secondary_electrical_conductivity = electrical_conductivity
    user_electrical_contact_conductance = 2.5e5 # as described in Cincotti et al (DOI: 10.1002/aic.11102)
    user_thermal_contact_conductance = 7 # also from Cincotti et al
    boundary = ssg_interface
  []
  [electric_contact_conductance]
    type = ElectrostaticContactCondition
    variable = potential_stainless_steel
    neighbor_var = potential_graphite
    primary_conductivity = electrical_conductivity
    secondary_conductivity = electrical_conductivity
    boundary = ssg_interface
    user_electrical_contact_conductance = 2.5e5 # as described in Cincotti et al (DOI: 10.1002/aic.11102)
  []

  [thermal_contact_conductance_calculated]
    type = ThermalContactCondition
    variable = temperature_stainless_steel
    neighbor_var = temperature_graphite
    primary_potential = potential_stainless_steel
    secondary_potential = potential_graphite
    primary_thermal_conductivity = thermal_conductivity
    secondary_thermal_conductivity = thermal_conductivity
    primary_electrical_conductivity = electrical_conductivity
    secondary_electrical_conductivity = electrical_conductivity
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = 8.52842e10 # resulting in electrical contact conductance = ~1.4715e5, thermal contact conductance = ~3.44689e7
    boundary = ssg_interface
  []
[]

[Materials]
  active = 'heat_conductor_graphite rho_graphite sigma_graphite heat_conductor_stainless_steel rho_stainless_steel sigma_stainless_steel'

  #graphite
  [heat_conductor_graphite]
    type = ADHeatConductionMaterial
    thermal_conductivity = 630
    specific_heat = 60
    block = graphite
  []
  [rho_graphite]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = 1.75e3
    block = graphite
  []
  [sigma_graphite]
    type = ADElectricalConductivity
    temperature = temperature_graphite
    reference_temperature = 293.0
    reference_resistivity = 3.0e-3
    temperature_coefficient = 0 # makes conductivity constant
    block = graphite
  []

  #stainless_steel
  [heat_conductor_stainless_steel]
    type = ADHeatConductionMaterial
    thermal_conductivity = 17
    specific_heat = 502
    block = stainless_steel
  []
  [rho_stainless_steel]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = 8e3
    block = stainless_steel
  []
  [sigma_stainless_steel]
    type = ADElectricalConductivity
    temperature = temperature_stainless_steel
    reference_temperature = 293.0
    reference_resistivity = 7e-7
    temperature_coefficient = 0 # makes conductivity constant
    block = stainless_steel
  []

  # harmonic mean of graphite and stainless steel hardness
  [mean_hardness]
    type = GraphiteStainlessMeanHardness
  []
[]

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

[Executioner]
  type = Transient
  scheme = bdf2
  solve_type = NEWTON
  dt = 0.1
  end_time = 100
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
[]

[Outputs]
  exodus = true
  perf_graph = true
[]

(test/tests/interfacekernels/thermal_conductance/thermal_interface.i)

# Input file block to generate this mesh can be found in thermal_interface_mesh.i
[Mesh]
  file = thermal_interface_regular_mesh.e
  coord_type = RZ
[]

[Variables]
  [temperature_graphite]
    initial_condition = 300.0
    block = graphite
  []
  [temperature_stainless_steel]
    initial_condition = 300.0
    block = stainless_steel
  []
  [potential_graphite]
    block = graphite
  []
  [potential_stainless_steel]
    block = stainless_steel
  []
[]

[AuxVariables]
  [T_infinity]
    initial_condition = 300.0
  []
[]

[Kernels]
  [HeatDiff_graphite]
    type = ADHeatConduction
    variable = temperature_graphite
    block = graphite
  []
  [HeatTdot_graphite]
    type = ADHeatConductionTimeDerivative
    variable = temperature_graphite
    block = graphite
  []
  [HeatSource_graphite]
    type = ADJouleHeatingSource
    variable = temperature_graphite
    elec = potential_graphite
    electrical_conductivity = electrical_conductivity
    block = graphite
  []

  [HeatDiff_stainless_steel]
    type = ADHeatConduction
    variable = temperature_stainless_steel
    block = stainless_steel
  []
  [HeatTdot_stainless_steel]
    type = ADHeatConductionTimeDerivative
    variable = temperature_stainless_steel
    block = stainless_steel
  []
  [HeatSource_stainless_steel]
    type = ADJouleHeatingSource
    variable = temperature_stainless_steel
    elec = potential_stainless_steel
    electrical_conductivity = electrical_conductivity
    block = stainless_steel
  []

  [electric_graphite]
    type = ADMatDiffusion
    variable = potential_graphite
    diffusivity = electrical_conductivity
    block = graphite
  []
  [electric_stainless_steel]
    type = ADMatDiffusion
    variable = potential_stainless_steel
    diffusivity = electrical_conductivity
    block = stainless_steel
  []
[]

[BCs]
  [external_surface_stainless]
    type = ADCoupledSimpleRadiativeHeatFluxBC
    boundary = right_stainless_steel
    variable = temperature_stainless_steel
    T_infinity = T_infinity
    emissivity = 0.85
  []
  [external_surface_graphite]
    type = ADCoupledSimpleRadiativeHeatFluxBC
    boundary = right_graphite
    variable = temperature_graphite
    T_infinity = T_infinity
    emissivity = 0.85
  []
  [elec_top]
    type = ADDirichletBC
    variable = potential_stainless_steel
    boundary = top_die
    value = 0.68 # Better reflects Cincotti et al (DOI: 10.1002/aic.11102) Figure 19
  []
  [elec_bottom]
    type = ADDirichletBC
    variable = potential_stainless_steel
    boundary = bottom_die
    value = 0
  []
[]

[InterfaceKernels]
  active = 'thermal_contact_conductance electric_contact_conductance'

  [thermal_contact_conductance]
    type = ThermalContactCondition
    variable = temperature_stainless_steel
    neighbor_var = temperature_graphite
    primary_potential = potential_stainless_steel
    secondary_potential = potential_graphite
    primary_thermal_conductivity = thermal_conductivity
    secondary_thermal_conductivity = thermal_conductivity
    primary_electrical_conductivity = electrical_conductivity
    secondary_electrical_conductivity = electrical_conductivity
    user_electrical_contact_conductance = 2.5e5 # as described in Cincotti et al (DOI: 10.1002/aic.11102)
    user_thermal_contact_conductance = 7 # also from Cincotti et al
    boundary = ssg_interface
  []
  [electric_contact_conductance]
    type = ElectrostaticContactCondition
    variable = potential_stainless_steel
    neighbor_var = potential_graphite
    primary_conductivity = electrical_conductivity
    secondary_conductivity = electrical_conductivity
    boundary = ssg_interface
    user_electrical_contact_conductance = 2.5e5 # as described in Cincotti et al (DOI: 10.1002/aic.11102)
  []

  [thermal_contact_conductance_calculated]
    type = ThermalContactCondition
    variable = temperature_stainless_steel
    neighbor_var = temperature_graphite
    primary_potential = potential_stainless_steel
    secondary_potential = potential_graphite
    primary_thermal_conductivity = thermal_conductivity
    secondary_thermal_conductivity = thermal_conductivity
    primary_electrical_conductivity = electrical_conductivity
    secondary_electrical_conductivity = electrical_conductivity
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = 8.52842e10 # resulting in electrical contact conductance = ~1.4715e5, thermal contact conductance = ~3.44689e7
    boundary = ssg_interface
  []
[]

[Materials]
  active = 'heat_conductor_graphite rho_graphite sigma_graphite heat_conductor_stainless_steel rho_stainless_steel sigma_stainless_steel'

  #graphite
  [heat_conductor_graphite]
    type = ADHeatConductionMaterial
    thermal_conductivity = 630
    specific_heat = 60
    block = graphite
  []
  [rho_graphite]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = 1.75e3
    block = graphite
  []
  [sigma_graphite]
    type = ADElectricalConductivity
    temperature = temperature_graphite
    reference_temperature = 293.0
    reference_resistivity = 3.0e-3
    temperature_coefficient = 0 # makes conductivity constant
    block = graphite
  []

  #stainless_steel
  [heat_conductor_stainless_steel]
    type = ADHeatConductionMaterial
    thermal_conductivity = 17
    specific_heat = 502
    block = stainless_steel
  []
  [rho_stainless_steel]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = 8e3
    block = stainless_steel
  []
  [sigma_stainless_steel]
    type = ADElectricalConductivity
    temperature = temperature_stainless_steel
    reference_temperature = 293.0
    reference_resistivity = 7e-7
    temperature_coefficient = 0 # makes conductivity constant
    block = stainless_steel
  []

  # harmonic mean of graphite and stainless steel hardness
  [mean_hardness]
    type = GraphiteStainlessMeanHardness
  []
[]

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

[Executioner]
  type = Transient
  scheme = bdf2
  solve_type = NEWTON
  dt = 0.1
  end_time = 100
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
[]

[Outputs]
  exodus = true
  perf_graph = true
[]

(test/tests/example_testing/sample_three.i)

# Cincotti Sample III SPS System Electrothermal Example
# Configuration: Two stainless steel electrodes (with cooling channels) with two
#                large spacers, two small spacers, and one plunger in-between.
#                RZ, axisymmetric. Reproduces Figures 20(a) and 20(b) in Cincotti
#                reference (below)
#
# BCs:
#    Potential:
#       V (top electrode, top surface) --> Neumann condition specifiying potential
#                                          based on applied RMS Current (670 A)
#                                          and cross-sectional electrode area (see
#                                          elec_top BC below). Current turned off
#                                          at 1200 s.
#       V (bottom electrode, bottom surface) = 0 V
#       V (elsewhere) --> natural boundary conditions (no current external to circuit)
#    Temperature:
#       T (top electrode, top surface) = 293 K
#       T (bottom electrode, bottom surface) = 293 K
#       T (right side) --> simple radiative BC into black body at 293 K
#       T (water channel) --> simple convective BC into "water" with heat transfer coefficient parameter from Cincotti and temperature of 293 K
#       T (elsewhere) --> natural boundary conditions
# Interface Conditions:
#       V (SS-G upper and G-SS lower interface) --> ElectrostaticContactCondition
#       T (SS-G upper and G-SS lower interface) --> ThermalContactCondition
# Initial Conditions:
#       V = default (0 V)
#       T = 293 K
# Applied Mechanical Load: 3.5 kN
#
# Reference: Cincotti et al, DOI 10.1002/aic.11102

[Mesh]
  [import_mesh]
    type = FileMeshGenerator
    file = spsdie_table1model3_2d.e
  []
  coord_type = RZ
[]

[Variables]
  [temperature_graphite]
    initial_condition = 293.0
    block = graphite
  []
  [temperature_stainless_steel]
    initial_condition = 293.0
    block = stainless_steel
  []
  [potential_graphite]
    block = graphite
  []
  [potential_stainless_steel]
    block = stainless_steel
  []
[]

[AuxVariables]
  [electric_field_r]
    family = MONOMIAL
    order = CONSTANT
  []
  [electric_field_z]
    family = MONOMIAL
    order = CONSTANT
  []
  [current_density]
    family = NEDELEC_ONE
    order = FIRST
  []
  [T_infinity]
    initial_condition = 293.0
  []
  [heatflux_graphite_r]
    family = MONOMIAL
    order = CONSTANT
    block = graphite
  []
  [heatflux_graphite_z]
    family = MONOMIAL
    order = CONSTANT
    block = graphite
  []
  [heatflux_stainless_steel_r]
    family = MONOMIAL
    order = CONSTANT
    block = stainless_steel
  []
  [heatflux_stainless_steel_z]
    family = MONOMIAL
    order = CONSTANT
    block = stainless_steel
  []
[]

[Kernels]
  [HeatDiff_graphite]
    type = ADHeatConduction
    variable = temperature_graphite
    block = graphite
  []
  [HeatTdot_graphite]
    type = ADHeatConductionTimeDerivative
    variable = temperature_graphite
    specific_heat = heat_capacity
    block = graphite
  []
  [HeatSource_graphite]
    type = ADJouleHeatingSource
    variable = temperature_graphite
    elec = potential_graphite
    block = graphite
  []

  [HeatDiff_stainless_steel]
    type = ADHeatConduction
    variable = temperature_stainless_steel
    block = stainless_steel
  []
  [HeatTdot_stainless_steel]
    type = ADHeatConductionTimeDerivative
    variable = temperature_stainless_steel
    specific_heat = heat_capacity
    block = stainless_steel
  []
  [HeatSource_stainless_steel]
    type = ADJouleHeatingSource
    variable = temperature_stainless_steel
    elec = potential_stainless_steel
    block = stainless_steel
  []

  [electric_graphite]
    type = ADMatDiffusion
    variable = potential_graphite
    diffusivity = electrical_conductivity
    block = graphite
  []
  [electric_stainless_steel]
    type = ADMatDiffusion
    variable = potential_stainless_steel
    diffusivity = electrical_conductivity
    block = stainless_steel
  []
[]

[AuxKernels]
  [electrostatic_calculation_r_graphite]
    type = PotentialToFieldAux
    gradient_variable = potential_graphite
    variable = electric_field_r
    sign = negative
    component = x
    block = graphite
  []
  [electrostatic_calculation_z_graphite]
    type = PotentialToFieldAux
    gradient_variable = potential_graphite
    variable = electric_field_z
    sign = negative
    component = y
    block = graphite
  []
  [electrostatic_calculation_r_stainless_steel]
    type = PotentialToFieldAux
    gradient_variable = potential_stainless_steel
    variable = electric_field_r
    sign = negative
    component = x
    block = stainless_steel
  []
  [electrostatic_calculation_z_stainless_steel]
    type = PotentialToFieldAux
    gradient_variable = potential_stainless_steel
    variable = electric_field_z
    sign = negative
    component = y
    block = stainless_steel
  []
  [current_density_graphite]
    type = ADCurrentDensity
    variable = current_density
    potential = potential_graphite
    block = graphite
  []
  [current_density_stainless_steel]
    type = ADCurrentDensity
    variable = current_density
    potential = potential_stainless_steel
    block = stainless_steel
  []
  [heat_flux_graphite_r]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature_graphite
    variable = heatflux_graphite_r
    component = x
    block = graphite
  []
  [heat_flux_graphite_z]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature_graphite
    variable = heatflux_graphite_z
    component = y
    block = graphite
  []
  [heat_flux_stainless_r]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature_stainless_steel
    variable = heatflux_stainless_steel_r
    block = stainless_steel
    component = x
  []
  [heat_flux_stainless_z]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature_stainless_steel
    variable = heatflux_stainless_steel_z
    block = stainless_steel
    component = y
  []
[]

[BCs]
  [external_surface_stainless]
    type = ADCoupledSimpleRadiativeHeatFluxBC
    boundary = right_die_stainless_steel
    variable = temperature_stainless_steel
    T_infinity = T_infinity
    emissivity = 0.4
  []
  [external_surface_graphite]
    type = ADCoupledSimpleRadiativeHeatFluxBC
    boundary = right_die_graphite
    variable = temperature_graphite
    T_infinity = T_infinity
    emissivity = 0.85
  []
  [water_channel]
    type = CoupledConvectiveHeatFluxBC
    boundary = water_channel
    variable = temperature_stainless_steel
    T_infinity = 293
    htc = 4725
  []
  [temp_top]
    type = ADDirichletBC
    variable = temperature_stainless_steel
    boundary = top_die
    value = 293
  []
  [temp_bottom]
    type = ADDirichletBC
    variable = temperature_stainless_steel
    boundary = bottom_die
    value = 293
  []
  [elec_top]
    type = ADFunctionNeumannBC
    variable = potential_stainless_steel
    boundary = top_die
    function = 'if(t < 47, (670 / (pi * 0.00155))*((0.057531/2.6969)*t), if(t > 1200, 0, 670 / (pi * 0.00155)))' # RMS Current / Cross-sectional Area. Ramping for t < 47s approximately reflects Cincotti et al (DOI: 10.1002/aic.11102) Figure 20(b)
  []
  [elec_bottom]
    type = ADDirichletBC
    variable = potential_stainless_steel
    boundary = bottom_die
    value = 0
  []
[]

[InterfaceKernels]
  [electric_contact_conductance_ssg]
    type = ElectrostaticContactCondition
    variable = potential_stainless_steel
    neighbor_var = potential_graphite
    boundary = ssg_interface
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
  []
  [thermal_contact_conductance_calculated_ssg]
    type = ThermalContactCondition
    variable = temperature_stainless_steel
    neighbor_var = temperature_graphite
    primary_potential = potential_stainless_steel
    secondary_potential = potential_graphite
    primary_thermal_conductivity = thermal_conductivity
    secondary_thermal_conductivity = thermal_conductivity
    primary_electrical_conductivity = electrical_conductivity
    secondary_electrical_conductivity = electrical_conductivity
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
    boundary = ssg_interface
  []
  [electric_contact_conductance_gss]
    type = ElectrostaticContactCondition
    variable = potential_graphite
    neighbor_var = potential_stainless_steel
    boundary = gss_interface
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
  []
  [thermal_contact_conductance_calculated_gss]
    type = ThermalContactCondition
    variable = temperature_graphite
    neighbor_var = temperature_stainless_steel
    primary_potential = potential_graphite
    secondary_potential = potential_stainless_steel
    primary_thermal_conductivity = thermal_conductivity
    secondary_thermal_conductivity = thermal_conductivity
    primary_electrical_conductivity = electrical_conductivity
    secondary_electrical_conductivity = electrical_conductivity
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
    boundary = gss_interface
  []
[]

[Materials]
  #graphite
  [heat_conductor_graphite]
    type = ADGraphiteThermal
    temperature = temperature_graphite
    block = graphite
  []
  [rho_graphite]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = 1750
    block = graphite
  []
  [sigma_graphite]
    type = ADGraphiteElectricalConductivity
    temperature = temperature_graphite
    block = graphite
  []

  #stainless_steel
  [heat_conductor_stainless_steel]
    type = ADStainlessSteelThermal
    temperature = temperature_stainless_steel
    block = stainless_steel
  []
  [rho_stainless_steel]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = 8000
    block = stainless_steel
  []
  [sigma_stainless_steel]
    type = ADStainlessSteelElectricalConductivity
    temperature = temperature_stainless_steel
    block = stainless_steel
  []

  # harmonic mean of graphite and stainless steel hardness
  [mean_hardness]
    type = GraphiteStainlessMeanHardness
  []

  # Material property converter for DiffusionFluxAux object
  [converter]
    type = MaterialADConverter
    ad_props_in = thermal_conductivity
    reg_props_out = nonad_thermal_conductivity
  []
[]

[Functions]
  [mechanical_pressure_func]
    type = ParsedFunction
    symbol_names = 'radius force'
    symbol_values = '0.04 3.5' # 'm kN'
    expression = 'force * 1e3 / (pi * radius^2)' # (N / m^2)
  []
[]

# Tracking data locations specified in Cincotti paper
# [Postprocessors]
#   [./lower_big_spacer]
#     type = PointValue
#     variable = temperature_graphite
#     point = '0 0.28 0'
#   [../]
#   [./lower_small_spacer]
#     type = PointValue
#     variable = temperature_graphite
#     point = '0 0.315 0'
#   [../]
#   [./center_plunger]
#     type = PointValue
#     variable = temperature_graphite
#     point = '0 0.34 0'
#   [../]
#   [./potential_tracking]
#     type = PointValue
#     variable = potential_stainless_steel
#     point = '0.027 0.584 0'
#   [../]
# []

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

[Executioner]
  type = Transient
  scheme = bdf2
  solve_type = NEWTON
  dt = 1
  end_time = 3600
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
[]

[Outputs]
  perf_graph = true
  file_base = sample_three_out
  active = 'normal_out'
  [normal_out]
    type = Exodus
  []
  [testing_out]
    # Due to platform diffs in *really* small auxvariable values
    type = Exodus
    show = 'temperature_graphite temperature_stainless_steel potential_graphite potential_stainless_steel'
  []
[]

(examples/sps/multiapp/complete_geometry/electrothermal/engineering_scale_electrothermal.i)

# Approximate Dr. Sinter Geometry
# Configuration: Two stainless steel electrodes (with cooling channels) with two
#                large spacers, and one die-and-plungers assembly.
#                RZ, axisymmetric.
# BCs:
#    Potential:
#       V (top electrode, top surface) --> Neumann condition specifiying potential
#                                          based on applied RMS Current (980 A)
#                                          and cross-sectional electrode area (see
#                                          elec_top BC below). Current turned off
#                                          at 1200 s.
#       V (bottom electrode, bottom surface) = 0 V
#       V (elsewhere) --> natural boundary conditions (no current external to circuit)
#    Temperature:
#       T (top electrode, top surface) = 293 K
#       T (bottom electrode, bottom surface) = 293 K
#       T (vertical right side) --> simple radiative BC into black body at 293 K
#       T (horizonal surfaces) --> GapHeatTransfer
#       T (water channel) --> simple convective BC into "water" with heat transfer coefficient parameter from Cincotti and temperature of 293 K
#       T (elsewhere) --> natural boundary conditions
#    Mechanics (not added yet):
#       Displacement: --> pinned axial centerpoint along interior surface of powder block, outer die wall
#       Pressure: top and bottom surfaces of upper and lower rams, 20.7 MPa
#                 ramped quickly at simulation start, then held constant
# Interface Conditions:
#       V (SS-G upper and G-SS lower interface) --> ElectrostaticContactCondition
#       T (SS-G upper and G-SS lower interface) --> ThermalContactCondition
#       V (graphite-yttria and graphite-graphite) --> GapHeatTransfer
#       T (graphite-yttria and graphite-graphite) --> GapHeatTransfer
#       P (graphite-yttria and graphite-graphite) --> Mechanical Contact
# Initial Conditions:
#       V = default (0 V)
#       T = 293 K
# Applied Mechanical Load: needs to be updated to 6.5kN, was 3.5 kN
#
# Reference for graphite, stainless steel: Cincotti et al, DOI 10.1002/aic.11102
# Assorted references for yttria, listed as comments in input file

initial_temperature = 293 #roughly 600C where the pyrometer kicks in
#initial_porosity=0.36 #Maximum random jammed packing, Donev et al (2004) Science Magazine

[Mesh]
  file = ../../../meshes/drsinter_nodalcontact_2d.e
  coord_type = RZ
  construct_side_list_from_node_list = true
  patch_update_strategy = iteration
  patch_size = 20
  # ghosting_patch_size = 5*patch_size
[]

[Problem]
  type = ReferenceResidualProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
[]

[Variables]
  [temperature_stainless_steel]
    initial_condition = ${initial_temperature}
    block = stainless_steel
  []
  [temperature]
    initial_condition = ${initial_temperature}
    block = 'graphite_spacers lower_plunger upper_plunger die_wall powder_compact'
  []
  [potential_stainless_steel]
    block = stainless_steel
  []
  [electric_potential]
    block = 'graphite_spacers lower_plunger upper_plunger die_wall powder_compact'
  []
[]

[AuxVariables]
  [electric_field_x]
    family = MONOMIAL #prettier pictures with smoother values
    order = FIRST
  []
  [electric_field_y]
    family = MONOMIAL
    order = FIRST
  []
  [current_density]
    family = NEDELEC_ONE
    order = FIRST
  []
  [yttria_sigma_aeh]
    initial_condition = 2.0e-10 #in units eV/((nV)^2-s-nm)
    block = 'powder_compact'
  []
  [microapp_potential]
    #converted to microapp electronVolts units
    block = 'powder_compact'
  []
  [E_x]
    order = FIRST
    family = MONOMIAL
    block = 'powder_compact'
  []
  [E_y]
    order = FIRST
    family = MONOMIAL
    block = 'powder_compact'
  []
  [yttria_current_density_forBC_microapp]
    order = FIRST
    family = MONOMIAL
    block = 'powder_compact'
  []

  # [T_infinity]
  #   initial_condition = ${initial_temperature}
  # []
  [heatflux_graphite_x]
    family = MONOMIAL
    order = CONSTANT
    block = 'graphite_spacers lower_plunger upper_plunger die_wall'
  []
  [heatflux_graphite_y]
    family = MONOMIAL
    order = CONSTANT
    block = 'graphite_spacers lower_plunger upper_plunger die_wall'
  []
  [heatflux_stainless_steel_x]
    family = MONOMIAL
    order = CONSTANT
    block = stainless_steel
  []
  [heatflux_stainless_steel_y]
    family = MONOMIAL
    order = CONSTANT
    block = stainless_steel
  []
  [heatflux_yttria_x]
    family = MONOMIAL
    order = CONSTANT
    block = 'powder_compact'
  []
  [heatflux_yttria_y]
    family = MONOMIAL
    order = CONSTANT
    block = 'powder_compact'
  []
  [yttria_thermal_conductivity_aeh]
    initial_condition = 0.4
  []
  [yttria_heat_capacity_volume_avg]
    initial_condition = 447.281 #842.2 at 1600K probably? # at 1500K #568.73 at 1000K #447.281 # at 293K
  []
  [yttria_density_volume_avg]
    initial_condition = 3206.4 ##5010.0*(1-${initial_porosity}) #in kg/m^3
  []
  [heat_transfer_radiation]
  []
  # [thermal_conductivity]
  #   family = MONOMIAL
  #   order = FIRST
  # []
[]

[Kernels]
  [HeatDiff_graphite]
    type = ADHeatConduction
    variable = temperature
    thermal_conductivity = thermal_conductivity
    extra_vector_tags = 'ref'
    block = 'graphite_spacers lower_plunger upper_plunger die_wall'
  []
  [HeatTdot_graphite]
    type = ADHeatConductionTimeDerivative
    variable = temperature
    specific_heat = heat_capacity
    density_name = density
    extra_vector_tags = 'ref'
    block = 'graphite_spacers lower_plunger upper_plunger die_wall'
  []
  [JouleHeating_graphite]
    type = ADJouleHeatingSource
    variable = temperature
    elec = electric_potential
    electrical_conductivity = electrical_conductivity
    use_displaced_mesh = true
    extra_vector_tags = 'ref'
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [electric_graphite]
    type = ADMatDiffusion
    variable = electric_potential
    diffusivity = electrical_conductivity
    extra_vector_tags = 'ref'
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []

  [HeatDiff_stainless_steel]
    type = ADHeatConduction
    variable = temperature_stainless_steel
    thermal_conductivity = thermal_conductivity
    extra_vector_tags = 'ref'
    block = stainless_steel
  []
  [HeatTdot_stainless_steel]
    type = ADHeatConductionTimeDerivative
    variable = temperature_stainless_steel
    specific_heat = heat_capacity
    density_name = density
    extra_vector_tags = 'ref'
    block = stainless_steel
  []
  [HeatSource_stainless_steel]
    type = ADJouleHeatingSource
    variable = temperature_stainless_steel
    elec = potential_stainless_steel
    electrical_conductivity = electrical_conductivity
    use_displaced_mesh = true
    extra_vector_tags = 'ref'
    block = stainless_steel
  []
  [electric_stainless_steel]
    type = ADMatDiffusion
    variable = potential_stainless_steel
    diffusivity = electrical_conductivity
    extra_vector_tags = 'ref'
    block = stainless_steel
  []

  [HeatDiff_yttria]
    type = ADHeatConduction
    variable = temperature
    thermal_conductivity = thermal_conductivity #use parsed material property, hope it works
    extra_vector_tags = 'ref'
    block = powder_compact
  []
  [HeatTdot_yttria]
    type = ADHeatConductionTimeDerivative
    variable = temperature
    specific_heat = heat_capacity #use parsed material property
    density_name = density
    extra_vector_tags = 'ref'
    block = powder_compact
  []
  [JouleHeating_yttria]
    type = ADJouleHeatingSource
    variable = temperature
    elec = electric_potential
    electrical_conductivity = electrical_conductivity
    use_displaced_mesh = true
    extra_vector_tags = 'ref'
    block = powder_compact
  []
  [electric_yttria]
    type = ADMatDiffusion
    variable = electric_potential
    diffusivity = electrical_conductivity
    extra_vector_tags = 'ref'
    block = powder_compact
  []
[]

[AuxKernels]
  [electrostatic_calculation_x_graphite_yttria]
    type = PotentialToFieldAux
    gradient_variable = electric_potential
    variable = electric_field_x
    sign = negative
    component = x
    block = 'powder_compact graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [electrostatic_calculation_y_graphite_yttria]
    type = PotentialToFieldAux
    gradient_variable = electric_potential
    variable = electric_field_y
    sign = negative
    component = y
    block = 'powder_compact graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [electrostatic_calculation_x_stainless_steel]
    type = PotentialToFieldAux
    gradient_variable = potential_stainless_steel
    variable = electric_field_x
    sign = negative
    component = x
    block = stainless_steel
  []
  [electrostatic_calculation_y_stainless_steel]
    type = PotentialToFieldAux
    gradient_variable = potential_stainless_steel
    variable = electric_field_y
    sign = negative
    component = y
    block = stainless_steel
  []
  [current_density_graphite_yttria]
    type = ADCurrentDensity
    variable = current_density
    potential = electric_potential
    block = 'powder_compact graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [current_density_stainless_steel]
    type = ADCurrentDensity
    variable = current_density
    potential = potential_stainless_steel
    block = stainless_steel
  []
  [microapp_potential]
    type = ParsedAux
    variable = microapp_potential
    coupled_variables = electric_potential
    expression = 'electric_potential*1e9' #convert from V to nV
    block = 'powder_compact'
  []
  [E_x]
    type = VariableGradientComponent
    variable = E_x
    gradient_variable = electric_potential
    component = x
    block = 'powder_compact'
  []
  [E_y]
    type = VariableGradientComponent
    variable = E_y
    gradient_variable = electric_potential
    component = y
    block = 'powder_compact'
  []
  [yttria_current_density_forBC_microapp]
    type = ParsedAux
    variable = yttria_current_density_forBC_microapp
    coupled_variables = 'electrical_conductivity E_y'
    expression = '-1.0*electrical_conductivity*E_y'
    block = 'powder_compact'
  []

  [heat_flux_graphite_x]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature
    variable = heatflux_graphite_x
    component = x
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [heat_flux_graphite_y]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature
    variable = heatflux_graphite_y
    component = y
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [heat_flux_stainless_x]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature_stainless_steel
    variable = heatflux_stainless_steel_x
    block = stainless_steel
    component = x
  []
  [heat_flux_stainless_y]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature_stainless_steel
    variable = heatflux_stainless_steel_y
    block = stainless_steel
    component = y
  []
  [heat_flux_yttria_x]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature
    variable = heatflux_yttria_x
    component = x
    block = 'powder_compact'
  []
  [heat_flux_yttria_y]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature
    variable = heatflux_yttria_y
    component = y
    block = 'powder_compact'
  []

  [heat_transfer_radiation_graphite]
    type = ParsedAux
    variable = heat_transfer_radiation
    boundary = 'outer_radiative_spacers outer_die_wall radiative_upper_plunger radiative_lower_plunger'
    coupled_variables = 'temperature'
    constant_names = 'boltzmann epsilon temperature_farfield' #published emissivity for graphite is 0.85
    constant_expressions = '5.67e-8 0.85 293.0' #roughly room temperature, which is probably too cold
    expression = '-boltzmann*epsilon*(temperature^4-temperature_farfield^4)'
  []
  [heat_transfer_radiation_stainless_steel]
    type = ParsedAux
    variable = heat_transfer_radiation
    boundary = 'outer_radiative_stainless_steel'
    coupled_variables = 'temperature_stainless_steel'
    constant_names = 'boltzmann epsilon temperature_farfield' #published emissivity for graphite is 0.85
    constant_expressions = '5.67e-8 0.4 293.0' #roughly room temperature, which is probably too cold
    expression = '-boltzmann*epsilon*(temperature_stainless_steel^4-temperature_farfield^4)'
  []
  [thermal_conductivity_graphite]
    type = ADMaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [thermal_conductivity_steel]
    type = ADMaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = 'stainless_steel'
  []
  [thermal_conductivity_yttria]
    type = ADMaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = 'powder_compact'
  []
  [electrical_conductivity_graphite]
    type = ADMaterialRealAux
    variable = electrical_conductivity
    property = electrical_conductivity
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [electrical_conductivity_steel]
    type = ADMaterialRealAux
    variable = electrical_conductivity
    property = electrical_conductivity
    block = 'stainless_steel'
  []
  [electrical_conductivity_yttria]
    type = ADMaterialRealAux
    variable = electrical_conductivity
    property = electrical_conductivity
    block = 'powder_compact'
  []
[]

[BCs]
  [external_surface_stainless]
    type = CoupledVarNeumannBC
    boundary = 'outer_radiative_stainless_steel'
    variable = temperature_stainless_steel
    v = heat_transfer_radiation
  []
  [external_surface_graphite]
    type = CoupledVarNeumannBC
    boundary = 'outer_radiative_spacers outer_die_wall radiative_upper_plunger radiative_lower_plunger'
    variable = temperature
    v = heat_transfer_radiation
  []
  [water_channel]
    type = CoupledConvectiveHeatFluxBC
    boundary = water_channel
    variable = temperature_stainless_steel
    T_infinity = ${initial_temperature}
    htc = 4725
  []
  [temperature_ram_extremes]
    type = ADDirichletBC
    variable = temperature_stainless_steel
    boundary = 'top_upper_steel_ram bottom_lower_steel_ram'
    value = 293
  []
  [electric_potential_top]
    type = ADFunctionNeumannBC
    variable = potential_stainless_steel
    boundary = top_upper_steel_ram
    ## This will need to be updated to match the Dr. Sinter geometry
    function = 'if(t < 31, (980 / (pi * 0.00155))*((0.141301/4.3625)*t), if(t > 600, 0, 980 / (pi * 0.00155)))' # RMS Current / Cross-sectional Area. Ramping for t < 31s approximately reflects Cincotti et al (DOI: 10.1002/aic.11102) Figure 21(b)
  []
  [electric_potential_bottom]
    type = ADDirichletBC
    variable = potential_stainless_steel
    boundary = bottom_lower_steel_ram
    value = 0
  []
[]

[Functions]
  [graphite_thermal_conductivity_fcn]
    type = ParsedFunction
    expression = '-4.418e-12*t^4+2.904e-8*t^3-4.688e-5*t^2-0.0316*t+119.659'
  []
  # [yttria_thermal_conductivity_fcn] #from the multiapp
  #   type = ParsedFunction
  #   expression = '3214.46/(t-147.73)'
  # []
  [harmonic_mean_thermal_conductivity]
    type = ParsedFunction
    expression = '2*(-4.418e-12*t^4+2.904e-8*t^3-4.688e-5*t^2-0.0316*t+119.659)*(3214.46/(t-147.73))/((-4.418e-12*t^4+2.904e-8*t^3-4.688e-5*t^2-0.0316*t+119.659)+(3214.46/(t-147.73)))'
    # symbol_names = 'k_graphite k_yttria'
    # symbol_values = 'graphite_thermal_conductivity_fcn yttria_thermal_conductivity_fcn'
  []

  [graphite_electrical_conductivity_fcn]
    type = ParsedFunction
    expression = '1.0/(-2.705e-15*t^3+1.263e-11*t^2-1.836e-8*t+1.813e-5)'
  []
  # [electrical_conductivity_fcn]
  #   type = ParsedFunction
  #   # symbol_names = porosity
  #   # symbol_values = initial_porosity
  #   expression = '(1-0.62)*2.0025e4*exp(-1.61/8.617343e-5/t)'
  # []
  [harmonic_mean_electrical_conductivity]
    type = ParsedFunction
    expression = '2*(1.0/(-2.705e-15*t^3+1.263e-11*t^2-1.836e-8*t+1.813e-5))*((1)*2.0025e4*exp(-1.61/8.617343e-5/t))/((1.0/(-2.705e-15*t^3+1.263e-11*t^2-1.836e-8*t+1.813e-5))+((1)*2.0025e4*exp(-1.61/8.617343e-5/t)))'
    # symbol_names = 'k_graphite k_yttria'
    # symbol_values = 'graphite_thermal_conductivity_fcn yttria_thermal_conductivity_fcn'
  []
  [mechanical_pressure_func]
    type = ParsedFunction
    symbol_names = 'radius coolant_radius force'
    symbol_values = '0.04 7.071e-3 6.5' # 'm kN'
    expression = 'force * 1e3 / (pi * (radius^2 - coolant_radius^2))' # (N / m^2)
  []
[]

[InterfaceKernels]
  [electric_contact_conductance_ssg]
    type = ElectrostaticContactCondition
    variable = potential_stainless_steel
    neighbor_var = electric_potential
    boundary = 'upper_ram_spacer_interface lower_ram_spacer_interface'
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
  []
  [thermal_contact_conductance_calculated_ssg]
    type = ThermalContactCondition
    variable = temperature_stainless_steel
    neighbor_var = temperature
    primary_potential = potential_stainless_steel
    secondary_potential = electric_potential
    primary_thermal_conductivity = thermal_conductivity
    secondary_thermal_conductivity = thermal_conductivity
    primary_electrical_conductivity = electrical_conductivity
    secondary_electrical_conductivity = electrical_conductivity
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
    boundary = 'upper_ram_spacer_interface lower_ram_spacer_interface'
  []
  [electric_contact_conductance_gss]
    type = ElectrostaticContactCondition
    variable = electric_potential
    neighbor_var = potential_stainless_steel
    boundary = 'upper_ram_spacer_interface lower_ram_spacer_interface'
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
  []
  [thermal_contact_conductance_calculated_gss]
    type = ThermalContactCondition
    variable = temperature
    neighbor_var = temperature_stainless_steel
    primary_potential = electric_potential
    secondary_potential = potential_stainless_steel
    primary_thermal_conductivity = thermal_conductivity
    secondary_thermal_conductivity = thermal_conductivity
    primary_electrical_conductivity = electrical_conductivity
    secondary_electrical_conductivity = electrical_conductivity
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
    boundary = 'upper_ram_spacer_interface lower_ram_spacer_interface'
  []
[]

[ThermalContact]
  [upper_plunger_powder_electric]
    type = GapHeatTransfer
    primary = bottom_upper_plunger
    secondary = top_powder_compact
    variable = electric_potential
    quadrature = true
    emissivity_primary = 0.0 #not applicable for electric potential
    emissivity_secondary = 0.0 #not applicable for electric potential
    gap_geometry_type = PLATE
    # gap_conductivity = 1e-23 #rough harmonic mean
    gap_conductivity_function = 'harmonic_mean_electrical_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [powder_bottom_plunger_electric]
    type = GapHeatTransfer
    primary = top_lower_plunger
    secondary = bottom_powder_compact
    variable = electric_potential
    quadrature = true
    emissivity_primary = 0.0 #not applicable for electric potential
    emissivity_secondary = 0.0 #not applicable for electric potential
    gap_geometry_type = PLATE
    # gap_conductivity = 1e-23 #rough harmonic mean
    gap_conductivity_function = 'harmonic_mean_electrical_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [powder_die_electric]
    type = GapHeatTransfer
    primary = inner_die_wall
    secondary = outer_powder_compact
    variable = electric_potential
    quadrature = true
    emissivity_primary = 0.0 #0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.0 #0.85
    # gap_geometry_type = PLATE # Not for vertical surfaces
    # gap_conductivity = 1e-23 #rough harmonic mean
    gap_conductivity_function = 'harmonic_mean_electrical_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [combined_plungers_die_electric]
    type = GapHeatTransfer
    primary = inner_die_wall
    secondary = combined_outer_plungers
    variable = electric_potential
    quadrature = true
    emissivity_primary = 0.0 #not applicable for electric potential
    emissivity_secondary = 0.0 #not applicable for electric potential
    # gap_geometry_type = PLATE # Not for vertical surfaces
    # gap_conductivity = 8.5e4 #graphite, at 500K
    gap_conductivity_function = 'graphite_electrical_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []

##Thermal Contact between gapped graphite die components
  [upper_plunger_spacer_gap_thermal]
    type = GapHeatTransfer
    primary = spacer_facing_upper_plunger
    secondary = plunger_facing_upper_spacer
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    gap_geometry_type = PLATE
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [upper_diewall_spacer_gap_thermal]
    type = GapHeatTransfer
    primary = top_die_wall
    secondary = bottom_upper_spacer
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    gap_geometry_type = PLATE
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [lower_plunger_spacer_gap_thermal]
    type = GapHeatTransfer
    primary = spacer_facing_lower_plunger
    secondary = plunger_facing_lower_spacer
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    gap_geometry_type = PLATE
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [lower_diewall_spacer_gap_thermal]
    type = GapHeatTransfer
    primary = bottom_die_wall
    secondary = top_lower_spacer
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    gap_geometry_type = PLATE
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [upper_plunger_diewall_gap_thermal]
    type = GapHeatTransfer
    primary = inner_die_wall ### paired temperature doesn't show on inner die wall, but temperature profile looks reasonable
    secondary = die_wall_facing_upper_plunger
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    # gap_geometry_type = PLATE # Not for vertical surfaces
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [lower_plunger_diewall_gap_thermal]
    type = GapHeatTransfer
    primary = inner_die_wall ### paired temperature doesn't show on inner die wall, but temperature profile looks reasonable
    secondary = die_wall_facing_lower_plunger
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    # gap_geometry_type = PLATE # Not for vertical surfaces
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []

  ## Thermal Contact between touching components of powder and die
  [upper_plunger_powder_thermal]
    type = GapHeatTransfer
    primary = bottom_upper_plunger
    secondary = top_powder_compact
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.3 #estimated from McMahon and Wilder, High Temperature Spectral Emissivity of Yttrium, Samarium, Gadolinium, Ebrium and Lutetium Oxides (1963) Atomic Energy Commission, IS-578, Figure 12
    gap_geometry_type = PLATE
    gap_conductivity_function = 'harmonic_mean_thermal_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [powder_bottom_plunger_thermal]
    type = GapHeatTransfer
    primary = top_lower_plunger
    secondary = bottom_powder_compact #expect more heat transfer from the die to the powder
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.3 #estimated from McMahon and Wilder, High Temperature Spectral Emissivity of Yttrium, Samarium, Gadolinium, Ebrium and Lutetium Oxides (1963) Atomic Energy Commission, IS-578, Figure 12
    gap_geometry_type = PLATE
    gap_conductivity_function = 'harmonic_mean_thermal_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [powder_die_thermal]
    type = GapHeatTransfer
    primary = inner_die_wall
    secondary = outer_powder_compact
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.3 #estimated from McMahon and Wilder, High Temperature Spectral Emissivity of Yttrium, Samarium, Gadolinium, Ebrium and Lutetium Oxides (1963) Atomic Energy Commission, IS-578, Figure 12
    # gap_geometry_type = PLATE # Not for vertical surfaces
    gap_conductivity_function = 'harmonic_mean_thermal_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [combined_plungers_die_thermal]
    type = GapHeatTransfer
    primary = inner_die_wall
    secondary = combined_outer_plungers
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    # gap_geometry_type = PLATE # Not for vertical surfaces
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
[]

[Materials]
  ## graphite blocks
  [graphite_thermal]
    type = ADGraphiteThermal
    temperature = temperature
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [graphite_density]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = 1.750e3 #in kg/m^3 from Cincotti et al 2007, Table 2, doi:10.1002/aic
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [graphite_electrical_conductivity]
    type = ADGraphiteElectricalConductivity
    temperature = temperature
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []

  #stainless_steel
  [stainless_steel_thermal]
    type = ADStainlessSteelThermal
    temperature = temperature_stainless_steel
    block = stainless_steel
  []
  [stainless_steel_density]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = 8000
    block = stainless_steel
  []
  [stainless_steel_electrical_conductivity]
    type = ADStainlessSteelElectricalConductivity
    temperature = temperature_stainless_steel
    block = stainless_steel
  []

  # harmonic mean of graphite and stainless steel hardness
  [mean_hardness_graphite_stainless_steel]
    type = GraphiteStainlessMeanHardness
    block = 'graphite_spacers upper_plunger lower_plunger die_wall stainless_steel'
  []

  ## yttria powder compact
  [yttria_thermal_conductivity]
    type = ADParsedMaterial
    coupled_variables = 'temperature'
    expression = '3214.46 / (temperature - 147.73)' #in W/(m-K) #Given from Larry's curve fitting, data from Klein and Croft, JAP, v. 38, p. 1603 and UC report "For Computer Heat Conduction Calculations - A compilation of thermal properties data" by A.L. Edwards, UCRL-50589 (1969)
    # coupled_variables = 'thermal_conductivity_aeh'
    # expression = 'thermal_conductivity_aeh' #in W/(m-K) directly, for now
    property_name = 'thermal_conductivity'
    output_properties = thermal_conductivity
    outputs = 'exodus'
    block = powder_compact
  []
  [yttria_specific_heat_capacity]
    type = ADParsedMaterial
    property_name = heat_capacity
    coupled_variables = 'yttria_heat_capacity_volume_avg'
    expression = 'yttria_heat_capacity_volume_avg' #in J/(K-kg)
    # output_properties = yttria_specific_heat_capacity
    # outputs = 'exodus'
    block = powder_compact
  []
  [yttria_density]
    type = ADParsedMaterial
    property_name = 'density'
    coupled_variables = 'yttria_density_volume_avg'
    expression = 'yttria_density_volume_avg'
    # output_properties = yttria_density
    # outputs = 'exodus'
    block = powder_compact
  []
  [electrical_conductivity]
    type = ADParsedMaterial
    #   coupled_variables = 'yttria_sigma_aeh'
    #   expression = 'yttria_sigma_aeh*1.602e8' #converts to units of J/(V^2-m-s)
    property_name = 'electrical_conductivity'
    output_properties = electrical_conductivity
    outputs = 'exodus'
    block = powder_compact
    # type = ADDerivativeParsedMaterial
    # property_name = electrical_conductivity
    coupled_variables = 'temperature'
    constant_names = 'Q_elec  kB            prefactor_solid  initial_porosity'
    constant_expressions = '1.61    8.617343e-5        1.25e-4           0.38'
    expression = '(1-initial_porosity) * prefactor_solid * exp(-Q_elec/kB/temperature) * 1.602e8' # in eV/(nV^2 s nm) per chat with Larry, last term converts to units of J/(V^2-m-s)
  []

  # Material property converter for DiffusionFluxAux object
  [converter]
    type = MaterialADConverter
    ad_props_in = thermal_conductivity
    reg_props_out = nonad_thermal_conductivity
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  automatic_scaling = true
  line_search = 'none'
  # compute_scaling_once = false

  # force running options
  # petsc_options_iname = '-pc_type -snes_linesearch_type -pc_factor_shift_type -pc_factor_shift_amount'
  # petsc_options_value = 'lu       basic                 NONZERO               1e-15'

  # #mechanical contact options
  # petsc_options = '-snes_ksp_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = ' lu       superlu_dist'
  petsc_options = '-snes_converged_reason -ksp_converged_reason'

  nl_forced_its = 1
  nl_rel_tol = 2e-5 #1e-6 #2e-5 for with mechanics #was 1e-10, for temperature only
  nl_abs_tol = 2e-12 #was 1e-12
  nl_max_its = 20
  l_max_its = 50
  dtmin = 1.0e-4

  end_time = 900 #600 #900 #15 minutes, rule of thumb from Dennis is 10 minutes
  [Quadrature]
    order = FIFTH #required for thermal and mechanical node-face contact
    side_order = SEVENTH
  []
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 0.05
    optimal_iterations = 8
    iteration_window = 2
  []
[]

[Postprocessors]
  [temperature_pp]
    type = AverageNodalVariableValue
    variable = temperature
    block = 'graphite_spacers upper_plunger lower_plunger die_wall powder_compact'
  []
  [temperature_stainless_steel_pp]
    type = AverageNodalVariableValue
    variable = temperature_stainless_steel
    block = stainless_steel
  []
  [graphite_thermal_conductivity]
    type = ElementAverageValue
    variable = thermal_conductivity
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [graphite_electrical_conductivity]
    type = ElementAverageValue
    variable = electrical_conductivity
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [steel_thermal_conductivity]
    type = ElementAverageValue
    variable = thermal_conductivity
    block = stainless_steel
  []
  [steel_electrical_conductivity]
    type = ElementAverageValue
    variable = electrical_conductivity
    block = stainless_steel
  []
  [yttria_thermal_conductivity]
    type = ElementAverageValue
    variable = thermal_conductivity
    block = powder_compact
  []
  [yttria_electrical_conductivity]
    type = ElementAverageValue
    variable = electrical_conductivity
    block = powder_compact
  []
[]

[Outputs]
  csv = true
  exodus = true
  perf_graph = true
  # [ckpt]
  #   type =Checkpoint
  #   time_step_interval = 1
  #   num_files = 2
  # []
[]

(test/tests/example_testing/sample_four.i)

# Cincotti Sample IV SPS System Electrothermal Example
# Configuration: Two stainless steel electrodes (with cooling channels) with two
#                large spacers, two small spacers, and one die-and-plungers
#                ensemble in-between. RZ, axisymmetric. Reproduces Figures 21(a)
#                and 21(b) in Cincotti reference (below)
#
# BCs:
#    Potential:
#       V (top electrode, top surface) --> Neumann condition specifiying potential
#                                          based on applied RMS Current (980 A)
#                                          and cross-sectional electrode area (see
#                                          elec_top BC below). Current turned off
#                                          at 1200 s.
#       V (bottom electrode, bottom surface) = 0 V
#       V (elsewhere) --> natural boundary conditions (no current external to circuit)
#    Temperature:
#       T (top electrode, top surface) = 293 K
#       T (bottom electrode, bottom surface) = 293 K
#       T (right side) --> simple radiative BC into black body at 293 K
#       T (water channel) --> simple convective BC into "water" with heat transfer coefficient parameter from Cincotti and temperature of 293 K
#       T (elsewhere) --> natural boundary conditions
# Interface Conditions:
#       V (SS-G upper and G-SS lower interface) --> ElectrostaticContactCondition
#       T (SS-G upper and G-SS lower interface) --> ThermalContactCondition
# Initial Conditions:
#       V = default (0 V)
#       T = 293 K
# Applied Mechanical Load: 3.5 kN
#
# Reference: Cincotti et al, DOI 10.1002/aic.11102

[Mesh]
  [import_mesh]
    type = FileMeshGenerator
    file = spsdie_table1model4_2d.e
  []
  coord_type = RZ
[]

[Variables]
  [temperature_graphite]
    initial_condition = 293.0
    block = graphite
  []
  [temperature_stainless_steel]
    initial_condition = 293.0
    block = stainless_steel
  []
  [potential_graphite]
    block = graphite
  []
  [potential_stainless_steel]
    block = stainless_steel
  []
[]

[AuxVariables]
  [electric_field_r]
    family = MONOMIAL
    order = CONSTANT
  []
  [electric_field_z]
    family = MONOMIAL
    order = CONSTANT
  []
  [current_density]
    family = NEDELEC_ONE
    order = FIRST
  []
  [T_infinity]
    initial_condition = 293.0
  []
  [heatflux_graphite_r]
    family = MONOMIAL
    order = CONSTANT
    block = graphite
  []
  [heatflux_graphite_z]
    family = MONOMIAL
    order = CONSTANT
    block = graphite
  []
  [heatflux_stainless_steel_r]
    family = MONOMIAL
    order = CONSTANT
    block = stainless_steel
  []
  [heatflux_stainless_steel_z]
    family = MONOMIAL
    order = CONSTANT
    block = stainless_steel
  []
[]

[Kernels]
  [HeatDiff_graphite]
    type = ADHeatConduction
    variable = temperature_graphite
    block = graphite
  []
  [HeatTdot_graphite]
    type = ADHeatConductionTimeDerivative
    variable = temperature_graphite
    specific_heat = heat_capacity
    block = graphite
  []
  [HeatSource_graphite]
    type = ADJouleHeatingSource
    variable = temperature_graphite
    elec = potential_graphite
    block = graphite
  []

  [HeatDiff_stainless_steel]
    type = ADHeatConduction
    variable = temperature_stainless_steel
    block = stainless_steel
  []
  [HeatTdot_stainless_steel]
    type = ADHeatConductionTimeDerivative
    variable = temperature_stainless_steel
    specific_heat = heat_capacity
    block = stainless_steel
  []
  [HeatSource_stainless_steel]
    type = ADJouleHeatingSource
    variable = temperature_stainless_steel
    elec = potential_stainless_steel
    block = stainless_steel
  []

  [electric_graphite]
    type = ADMatDiffusion
    variable = potential_graphite
    diffusivity = electrical_conductivity
    block = graphite
  []
  [electric_stainless_steel]
    type = ADMatDiffusion
    variable = potential_stainless_steel
    diffusivity = electrical_conductivity
    block = stainless_steel
  []
[]

[AuxKernels]
  [electrostatic_calculation_r_graphite]
    type = PotentialToFieldAux
    gradient_variable = potential_graphite
    variable = electric_field_r
    sign = negative
    component = x
    block = graphite
  []
  [electrostatic_calculation_z_graphite]
    type = PotentialToFieldAux
    gradient_variable = potential_graphite
    variable = electric_field_z
    sign = negative
    component = y
    block = graphite
  []
  [electrostatic_calculation_r_stainless_steel]
    type = PotentialToFieldAux
    gradient_variable = potential_stainless_steel
    variable = electric_field_r
    sign = negative
    component = x
    block = stainless_steel
  []
  [electrostatic_calculation_z_stainless_steel]
    type = PotentialToFieldAux
    gradient_variable = potential_stainless_steel
    variable = electric_field_z
    sign = negative
    component = y
    block = stainless_steel
  []
  [current_density_graphite]
    type = ADCurrentDensity
    variable = current_density
    potential = potential_graphite
    block = graphite
  []
  [current_density_stainless_steel]
    type = ADCurrentDensity
    variable = current_density
    potential = potential_stainless_steel
    block = stainless_steel
  []
  [heat_flux_graphite_r]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature_graphite
    variable = heatflux_graphite_r
    component = x
    block = graphite
  []
  [heat_flux_graphite_z]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature_graphite
    variable = heatflux_graphite_z
    component = y
    block = graphite
  []
  [heat_flux_stainless_r]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature_stainless_steel
    variable = heatflux_stainless_steel_r
    block = stainless_steel
    component = x
  []
  [heat_flux_stainless_z]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature_stainless_steel
    variable = heatflux_stainless_steel_z
    block = stainless_steel
    component = y
  []
[]

[BCs]
  [external_surface_stainless]
    type = ADCoupledSimpleRadiativeHeatFluxBC
    boundary = right_die_stainless_steel
    variable = temperature_stainless_steel
    T_infinity = T_infinity
    emissivity = 0.4
  []
  [external_surface_graphite]
    type = ADCoupledSimpleRadiativeHeatFluxBC
    boundary = right_die_graphite
    variable = temperature_graphite
    T_infinity = T_infinity
    emissivity = 0.85
  []
  [water_channel]
    type = CoupledConvectiveHeatFluxBC
    boundary = water_channel
    variable = temperature_stainless_steel
    T_infinity = 293
    htc = 4725
  []
  [temp_top]
    type = ADDirichletBC
    variable = temperature_stainless_steel
    boundary = top_die
    value = 293
  []
  [temp_bottom]
    type = ADDirichletBC
    variable = temperature_stainless_steel
    boundary = bottom_die
    value = 293
  []
  [elec_top]
    type = ADFunctionNeumannBC
    variable = potential_stainless_steel
    boundary = top_die
    function = 'if(t < 31, (980 / (pi * 0.00155))*((0.141301/4.3625)*t), if(t > 600, 0, 980 / (pi * 0.00155)))' # RMS Current / Cross-sectional Area. Ramping for t < 31s approximately reflects Cincotti et al (DOI: 10.1002/aic.11102) Figure 21(b)
  []
  [elec_bottom]
    type = ADDirichletBC
    variable = potential_stainless_steel
    boundary = bottom_die
    value = 0
  []
[]

[InterfaceKernels]
  [electric_contact_conductance_ssg]
    type = ElectrostaticContactCondition
    variable = potential_stainless_steel
    neighbor_var = potential_graphite
    boundary = ssg_interface
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
  []
  [thermal_contact_conductance_calculated_ssg]
    type = ThermalContactCondition
    variable = temperature_stainless_steel
    neighbor_var = temperature_graphite
    primary_potential = potential_stainless_steel
    secondary_potential = potential_graphite
    primary_thermal_conductivity = thermal_conductivity
    secondary_thermal_conductivity = thermal_conductivity
    primary_electrical_conductivity = electrical_conductivity
    secondary_electrical_conductivity = electrical_conductivity
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
    boundary = ssg_interface
  []
  [electric_contact_conductance_gss]
    type = ElectrostaticContactCondition
    variable = potential_graphite
    neighbor_var = potential_stainless_steel
    boundary = gss_interface
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
  []
  [thermal_contact_conductance_calculated_gss]
    type = ThermalContactCondition
    variable = temperature_graphite
    neighbor_var = temperature_stainless_steel
    primary_potential = potential_graphite
    secondary_potential = potential_stainless_steel
    primary_thermal_conductivity = thermal_conductivity
    secondary_thermal_conductivity = thermal_conductivity
    primary_electrical_conductivity = electrical_conductivity
    secondary_electrical_conductivity = electrical_conductivity
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
    boundary = gss_interface
  []
[]

[Materials]
  #graphite
  [heat_conductor_graphite]
    type = ADGraphiteThermal
    temperature = temperature_graphite
    block = graphite
  []
  [rho_graphite]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = 1750
    block = graphite
  []
  [sigma_graphite]
    type = ADGraphiteElectricalConductivity
    temperature = temperature_graphite
    block = graphite
  []

  #stainless_steel
  [heat_conductor_stainless_steel]
    type = ADStainlessSteelThermal
    temperature = temperature_stainless_steel
    block = stainless_steel
  []
  [rho_stainless_steel]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = 8000
    block = stainless_steel
  []
  [sigma_stainless_steel]
    type = ADStainlessSteelElectricalConductivity
    temperature = temperature_stainless_steel
    block = stainless_steel
  []

  # harmonic mean of graphite and stainless steel hardness
  [mean_hardness]
    type = GraphiteStainlessMeanHardness
  []

  # Material property converter for DiffusionFluxAux object
  [converter]
    type = MaterialADConverter
    ad_props_in = thermal_conductivity
    reg_props_out = nonad_thermal_conductivity
  []
[]

[Functions]
  [mechanical_pressure_func]
    type = ParsedFunction
    symbol_names = 'radius force'
    symbol_values = '0.04 3.5' # 'm kN'
    expression = 'force * 1e3 / (pi * radius^2)' # (N / m^2)
  []
[]

# Tracking data locations specified in Cincotti paper
# [Postprocessors]
#   [./lower_small_spacer]
#     type = PointValue
#     variable = temperature_graphite
#     point = '0 0.315 0'
#   [../]
#   [./inside_die]
#     type = PointValue
#     variable = temperature_graphite
#     point = '0.011 0.3525 0'
#   [../]
#   [./potential_tracking]
#     type = PointValue
#     variable = potential_stainless_steel
#     point = '0.027 0.609 0'
#   [../]
# []

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

[Executioner]
  type = Transient
  scheme = bdf2
  solve_type = NEWTON
  dt = 1
  end_time = 1800
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
[]

[Outputs]
  perf_graph = true
  file_base = sample_four_out
  active = 'normal_out'
  [normal_out]
    type = Exodus
  []
  [testing_out]
    # Due to platform diffs in *really* small auxvariable values
    type = Exodus
    show = 'temperature_graphite temperature_stainless_steel potential_graphite potential_stainless_steel'
  []
[]

(test/tests/example_testing/sample_one.i)

# Cincotti Sample I SPS System Electrothermal Example
# Configuration: Two stainless steel electrodes (with cooling channels) with two
#                large spacers in-between. RZ, axisymmetric. Reproduces Figures
#                18(a), 18(b), and 19 in Cincotti reference (below)
#
# BCs:
#    Potential:
#       V (top electrode, top surface) --> Neumann condition specifiying potential
#                                          based on applied RMS Current (1200 A)
#                                          and cross-sectional electrode area (see
#                                          elec_top BC below). Current turned off
#                                          at 1200 s.
#       V (bottom electrode, bottom surface) = 0 V
#       V (elsewhere) --> natural boundary conditions (no current external to circuit)
#    Temperature:
#       T (top electrode, top surface) = 293 K
#       T (bottom electrode, bottom surface) = 293 K
#       T (right side) --> simple radiative BC into black body at 293 K
#       T (water channel) --> simple convective BC into "water" with heat transfer
#                             coefficient parameter from Cincotti and temperature of 293 K
#       T (elsewhere) --> natural boundary conditions
# Interface Conditions:
#       V (SS-G upper and G-SS lower interface) --> ElectrostaticContactCondition
#       T (SS-G upper and G-SS lower interface) --> ThermalContactCondition
# Initial Conditions:
#       V = default (0 V)
#       T = 293 K
# Applied Mechanical Load: 50 kN
#
# Reference: Cincotti et al, DOI 10.1002/aic.11102

[Mesh]
  [import_mesh]
    type = FileMeshGenerator
    file = spsdie_table1model1_2d.e
  []
  coord_type = RZ
[]

[Variables]
  [temperature_graphite]
    initial_condition = 293.0
    block = graphite
  []
  [temperature_stainless_steel]
    initial_condition = 293.0
    block = stainless_steel
  []
  [potential_graphite]
    block = graphite
  []
  [potential_stainless_steel]
    block = stainless_steel
  []
[]

[AuxVariables]
  [electric_field_r]
    family = MONOMIAL
    order = CONSTANT
  []
  [electric_field_z]
    family = MONOMIAL
    order = CONSTANT
  []
  [current_density]
    family = NEDELEC_ONE
    order = FIRST
  []
  [T_infinity]
    initial_condition = 293.0
  []
  [heatflux_graphite_r]
    family = MONOMIAL
    order = CONSTANT
    block = graphite
  []
  [heatflux_graphite_z]
    family = MONOMIAL
    order = CONSTANT
    block = graphite
  []
  [heatflux_stainless_steel_r]
    family = MONOMIAL
    order = CONSTANT
    block = stainless_steel
  []
  [heatflux_stainless_steel_z]
    family = MONOMIAL
    order = CONSTANT
    block = stainless_steel
  []
[]

[Kernels]
  [HeatDiff_graphite]
    type = ADHeatConduction
    variable = temperature_graphite
    block = graphite
  []
  [HeatTdot_graphite]
    type = ADHeatConductionTimeDerivative
    variable = temperature_graphite
    specific_heat = heat_capacity
    block = graphite
  []
  [HeatSource_graphite]
    type = ADJouleHeatingSource
    variable = temperature_graphite
    elec = potential_graphite
    block = graphite
  []

  [HeatDiff_stainless_steel]
    type = ADHeatConduction
    variable = temperature_stainless_steel
    block = stainless_steel
  []
  [HeatTdot_stainless_steel]
    type = ADHeatConductionTimeDerivative
    variable = temperature_stainless_steel
    specific_heat = heat_capacity
    block = stainless_steel
  []
  [HeatSource_stainless_steel]
    type = ADJouleHeatingSource
    variable = temperature_stainless_steel
    elec = potential_stainless_steel
    block = stainless_steel
  []

  [electric_graphite]
    type = ADMatDiffusion
    variable = potential_graphite
    diffusivity = electrical_conductivity
    block = graphite
  []
  [electric_stainless_steel]
    type = ADMatDiffusion
    variable = potential_stainless_steel
    diffusivity = electrical_conductivity
    block = stainless_steel
  []
[]

[AuxKernels]
  [electrostatic_calculation_r_graphite]
    type = PotentialToFieldAux
    gradient_variable = potential_graphite
    variable = electric_field_r
    sign = negative
    component = x
    block = graphite
  []
  [electrostatic_calculation_z_graphite]
    type = PotentialToFieldAux
    gradient_variable = potential_graphite
    variable = electric_field_z
    sign = negative
    component = y
    block = graphite
  []
  [electrostatic_calculation_r_stainless_steel]
    type = PotentialToFieldAux
    gradient_variable = potential_stainless_steel
    variable = electric_field_r
    sign = negative
    component = x
    block = stainless_steel
  []
  [electrostatic_calculation_z_stainless_steel]
    type = PotentialToFieldAux
    gradient_variable = potential_stainless_steel
    variable = electric_field_z
    sign = negative
    component = y
    block = stainless_steel
  []
  [current_density_graphite]
    type = ADCurrentDensity
    variable = current_density
    potential = potential_graphite
    block = graphite
  []
  [current_density_stainless_steel]
    type = ADCurrentDensity
    variable = current_density
    potential = potential_stainless_steel
    block = stainless_steel
  []
  [heat_flux_graphite_r]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature_graphite
    variable = heatflux_graphite_r
    component = x
    block = graphite
  []
  [heat_flux_graphite_z]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature_graphite
    variable = heatflux_graphite_z
    component = y
    block = graphite
  []
  [heat_flux_stainless_r]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature_stainless_steel
    variable = heatflux_stainless_steel_r
    block = stainless_steel
    component = x
  []
  [heat_flux_stainless_z]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature_stainless_steel
    variable = heatflux_stainless_steel_z
    block = stainless_steel
    component = y
  []
[]

[BCs]
  [external_surface_stainless]
    type = ADCoupledSimpleRadiativeHeatFluxBC
    boundary = right_die_stainless_steel
    variable = temperature_stainless_steel
    T_infinity = T_infinity
    emissivity = 0.4
  []
  [external_surface_graphite]
    type = ADCoupledSimpleRadiativeHeatFluxBC
    boundary = right_die_graphite
    variable = temperature_graphite
    T_infinity = T_infinity
    emissivity = 0.85
  []
  [water_channel]
    type = CoupledConvectiveHeatFluxBC
    boundary = water_channel
    variable = temperature_stainless_steel
    T_infinity = 293
    htc = 4725
  []
  [temp_top]
    type = ADDirichletBC
    variable = temperature_stainless_steel
    boundary = top_die
    value = 293
  []
  [temp_bottom]
    type = ADDirichletBC
    variable = temperature_stainless_steel
    boundary = bottom_die
    value = 293
  []
  [elec_top]
    type = ADFunctionNeumannBC
    variable = potential_stainless_steel
    boundary = top_die
    function = 'if(t < 43, (1200 / (pi * 0.00155))*((0.017703/0.75549)*t), if(t > 1200, 0, 1200 / (pi * 0.00155)))' # RMS Current / Cross-sectional Area. Ramping for t < 43s approximately reflects Cincotti et al (DOI: 10.1002/aic.11102) Figure 18(b)
  []
  [elec_bottom]
    type = ADDirichletBC
    variable = potential_stainless_steel
    boundary = bottom_die
    value = 0
  []
[]

[InterfaceKernels]
  [electric_contact_conductance_ssg]
    type = ElectrostaticContactCondition
    variable = potential_stainless_steel
    neighbor_var = potential_graphite
    boundary = ssg_interface
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
  []
  [thermal_contact_conductance_calculated_ssg]
    type = ThermalContactCondition
    variable = temperature_stainless_steel
    neighbor_var = temperature_graphite
    primary_potential = potential_stainless_steel
    secondary_potential = potential_graphite
    primary_thermal_conductivity = thermal_conductivity
    secondary_thermal_conductivity = thermal_conductivity
    primary_electrical_conductivity = electrical_conductivity
    secondary_electrical_conductivity = electrical_conductivity
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
    boundary = ssg_interface
  []
  [electric_contact_conductance_gss]
    type = ElectrostaticContactCondition
    variable = potential_graphite
    neighbor_var = potential_stainless_steel
    boundary = gss_interface
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
  []
  [thermal_contact_conductance_calculated_gss]
    type = ThermalContactCondition
    variable = temperature_graphite
    neighbor_var = temperature_stainless_steel
    primary_potential = potential_graphite
    secondary_potential = potential_stainless_steel
    primary_thermal_conductivity = thermal_conductivity
    secondary_thermal_conductivity = thermal_conductivity
    primary_electrical_conductivity = electrical_conductivity
    secondary_electrical_conductivity = electrical_conductivity
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
    boundary = gss_interface
  []
[]

[Materials]
  #graphite
  [heat_conductor_graphite]
    type = ADGraphiteThermal
    temperature = temperature_graphite
    block = graphite
  []
  [rho_graphite]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = 1750
    block = graphite
  []
  [sigma_graphite]
    type = ADGraphiteElectricalConductivity
    temperature = temperature_graphite
    block = graphite
  []

  #stainless_steel
  [heat_conductor_stainless_steel]
    type = ADStainlessSteelThermal
    temperature = temperature_stainless_steel
    block = stainless_steel
  []
  [rho_stainless_steel]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = 8000
    block = stainless_steel
  []
  [sigma_stainless_steel]
    type = ADStainlessSteelElectricalConductivity
    temperature = temperature_stainless_steel
    block = stainless_steel
  []

  # harmonic mean of graphite and stainless steel hardness
  [mean_hardness]
    type = GraphiteStainlessMeanHardness
  []

  # Material property converter for DiffusionFluxAux object
  [converter]
    type = MaterialADConverter
    ad_props_in = thermal_conductivity
    reg_props_out = nonad_thermal_conductivity
  []
[]

[Functions]
  [mechanical_pressure_func]
    type = ParsedFunction
    symbol_names = 'radius force'
    symbol_values = '0.04 50' # 'm kN'
    expression = 'force * 1e3 / (pi * radius^2)' # (N / m^2)
  []
[]

# Tracking data locations specified in Cincotti paper
# [Postprocessors]
#   [./temp_tracking]
#     type = PointValue
#     variable = temperature_graphite
#     point = '0 0.28 0'
#   [../]
#   [./potential_tracking]
#     type = PointValue
#     variable = potential_stainless_steel
#     point = '0.027 0.504 0'
#   [../]
# []

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

[Executioner]
  type = Transient
  scheme = bdf2
  solve_type = NEWTON
  dt = 1
  end_time = 3600
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
[]

[Outputs]
  perf_graph = true
  file_base = sample_one_out
  active = 'normal_out'
  [normal_out]
    type = Exodus
  []
  [testing_out]
    # Due to platform diffs in *really* small auxvariable values
    type = Exodus
    show = 'temperature_graphite temperature_stainless_steel potential_graphite potential_stainless_steel'
  []
[]

(examples/sps/multiapp/electrothermal_with_phase_field/twoway_lots_of_particles_prototype/engineering_scale_electrothermal_twoway_lots_prototype.i)

# Approximate Dr. Sinter Geometry
# Configuration: Two stainless steel electrodes (with cooling channels) with two
#                large spacers, and one die-and-plungers assembly.
#                RZ, axisymmetric.
#
# Two-way coupling (thermal and electrical conductivities and integrated heat
#    source) of micro-scale app with lots of particles in one location.
#
# BCs:
#    Potential:
#       V (top electrode, top surface) --> Neumann condition specifiying potential
#                                          based on applied RMS Current (980 A)
#                                          and cross-sectional electrode area (see
#                                          elec_top BC below). Current turned off
#                                          at 1200 s.
#       V (bottom electrode, bottom surface) = 0 V
#       V (elsewhere) --> natural boundary conditions (no current external to circuit)
#    Temperature:
#       T (top electrode, top surface) = 293 K
#       T (bottom electrode, bottom surface) = 293 K
#       T (vertical right side) --> simple radiative BC into black body at 293 K
#       T (horizonal surfaces) --> GapHeatTransfer
#       T (water channel) --> simple convective BC into "water" with heat transfer coefficient parameter from Cincotti and temperature of 293 K
#       T (elsewhere) --> natural boundary conditions
#    Mechanics (not added yet):
#       Displacement: --> pinned axial centerpoint along interior surface of powder block, outer die wall
#       Pressure: top and bottom surfaces of upper and lower rams, 20.7 MPa
#                 ramped quickly at simulation start, then held constant
# Interface Conditions:
#       V (SS-G upper and G-SS lower interface) --> ElectrostaticContactCondition
#       T (SS-G upper and G-SS lower interface) --> ThermalContactCondition
#       V (graphite-yttria and graphite-graphite) --> GapHeatTransfer
#       T (graphite-yttria and graphite-graphite) --> GapHeatTransfer
#       P (graphite-yttria and graphite-graphite) --> Mechanical Contact
# Initial Conditions:
#       V = default (0 V)
#       T = 873 K
# Applied Mechanical Load: needs to be updated to 6.5kN, was 3.5 kN
# Location of micro-app: (x, y, z) = (0.00125, 0.034, 0)
#
# Reference for graphite, stainless steel: Cincotti et al, DOI 10.1002/aic.11102
# Assorted references for yttria, listed as comments in input file

initial_temperature = 873 #roughly 600C where the pyrometer kicks in
#initial_porosity=0.36 #Maximum random jammed packing, Donev et al (2004) Science Magazine

[Mesh]
  file = ../../../meshes/drsinter_nodalcontact_2d.e
  coord_type = RZ
  construct_side_list_from_node_list = true
  patch_update_strategy = iteration
  patch_size = 20
  # ghosting_patch_size = 5*patch_size
[]

[Problem]
  type = ReferenceResidualProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
[]

[Variables]
  [temperature_stainless_steel]
    initial_condition = 300 #units of K
    block = stainless_steel
  []
  [temperature]
    initial_condition = ${initial_temperature}
    block = 'graphite_spacers lower_plunger upper_plunger die_wall powder_compact'
  []
  [potential_stainless_steel]
    block = stainless_steel
  []
  [electric_potential]
    block = 'graphite_spacers lower_plunger upper_plunger die_wall powder_compact'
  []
[]

[AuxVariables]
  [electric_field_x]
    family = MONOMIAL #prettier pictures with smoother values
    order = FIRST
  []
  [electric_field_y]
    family = MONOMIAL
    order = FIRST
  []
  [current_density]
    family = NEDELEC_ONE
    order = FIRST
  []
  [yttria_sigma_aeh]
    initial_condition = 2.0e-10 #in units eV/((nV)^2-s-nm)
    block = 'powder_compact'
  []
  [microapp_potential]
    #converted to microapp electronVolts units
    block = 'powder_compact'
  []
  [E_x]
    order = FIRST
    family = MONOMIAL
    block = 'powder_compact'
  []
  [E_y]
    order = FIRST
    family = MONOMIAL
    block = 'powder_compact'
  []
  [yttria_current_density_forBC_microapp]
    order = FIRST
    family = MONOMIAL
    block = 'powder_compact'
  []

  [Q_from_sub]
    #this will be in eV/m/s, will need unit conversion to J/m^3/s based on phase-field domain size
    order = FIRST
    family = LAGRANGE
  []

  # [T_infinity]
  #   initial_condition = ${initial_temperature}
  # []
  [heatflux_graphite_x]
    family = MONOMIAL
    order = CONSTANT
    block = 'graphite_spacers lower_plunger upper_plunger die_wall'
  []
  [heatflux_graphite_y]
    family = MONOMIAL
    order = CONSTANT
    block = 'graphite_spacers lower_plunger upper_plunger die_wall'
  []
  [heatflux_stainless_steel_x]
    family = MONOMIAL
    order = CONSTANT
    block = stainless_steel
  []
  [heatflux_stainless_steel_y]
    family = MONOMIAL
    order = CONSTANT
    block = stainless_steel
  []
  [heatflux_yttria_x]
    family = MONOMIAL
    order = CONSTANT
    block = 'powder_compact'
  []
  [heatflux_yttria_y]
    family = MONOMIAL
    order = CONSTANT
    block = 'powder_compact'
  []
  [yttria_thermal_conductivity_aeh]
    initial_condition = 0.4
  []
  [yttria_heat_capacity_volume_avg]
    initial_condition = 447.281 #842.2 at 1600K probably? # at 1500K #568.73 at 1000K #447.281 # at 293K
  []
  [yttria_density_volume_avg]
    initial_condition = 3206.4 ##5010.0*(1-${initial_porosity}) #in kg/m^3
  []
  [heat_transfer_radiation]
  []
  # [thermal_conductivity]
  #   family = MONOMIAL
  #   order = FIRST
  # []
[]

[Kernels]
  [HeatDiff_graphite]
    type = ADHeatConduction
    variable = temperature
    thermal_conductivity = thermal_conductivity
    extra_vector_tags = 'ref'
    block = 'graphite_spacers lower_plunger upper_plunger die_wall'
  []
  [HeatTdot_graphite]
    type = ADHeatConductionTimeDerivative
    variable = temperature
    specific_heat = heat_capacity
    density_name = density
    extra_vector_tags = 'ref'
    block = 'graphite_spacers lower_plunger upper_plunger die_wall'
  []
  [JouleHeating_graphite]
    type = ADJouleHeatingSource
    variable = temperature
    elec = electric_potential
    electrical_conductivity = electrical_conductivity
    use_displaced_mesh = true
    extra_vector_tags = 'ref'
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [electric_graphite]
    type = ADMatDiffusion
    variable = electric_potential
    diffusivity = electrical_conductivity
    extra_vector_tags = 'ref'
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []

  [HeatDiff_stainless_steel]
    type = ADHeatConduction
    variable = temperature_stainless_steel
    thermal_conductivity = thermal_conductivity
    extra_vector_tags = 'ref'
    block = stainless_steel
  []
  [HeatTdot_stainless_steel]
    type = ADHeatConductionTimeDerivative
    variable = temperature_stainless_steel
    specific_heat = heat_capacity
    density_name = density
    extra_vector_tags = 'ref'
    block = stainless_steel
  []
  [HeatSource_stainless_steel]
    type = ADJouleHeatingSource
    variable = temperature_stainless_steel
    elec = potential_stainless_steel
    electrical_conductivity = electrical_conductivity
    use_displaced_mesh = true
    extra_vector_tags = 'ref'
    block = stainless_steel
  []
  [electric_stainless_steel]
    type = ADMatDiffusion
    variable = potential_stainless_steel
    diffusivity = electrical_conductivity
    extra_vector_tags = 'ref'
    block = stainless_steel
  []

  [HeatDiff_yttria]
    type = ADHeatConduction
    variable = temperature
    thermal_conductivity = thermal_conductivity #use parsed material property, hope it works
    extra_vector_tags = 'ref'
    block = powder_compact
  []
  [HeatTdot_yttria]
    type = ADHeatConductionTimeDerivative
    variable = temperature
    specific_heat = heat_capacity #use parsed material property
    density_name = density
    extra_vector_tags = 'ref'
    block = powder_compact
  []
  [JouleHeating_yttria]
    type = ADJouleHeatingSource
    variable = temperature
    elec = electric_potential
    electrical_conductivity = electrical_conductivity
    use_displaced_mesh = true
    extra_vector_tags = 'ref'
    block = powder_compact
  []
  [electric_yttria]
    type = ADMatDiffusion
    variable = electric_potential
    diffusivity = electrical_conductivity
    extra_vector_tags = 'ref'
    block = powder_compact
  []
[]

[AuxKernels]
  [electrostatic_calculation_x_graphite_yttria]
    type = PotentialToFieldAux
    gradient_variable = electric_potential
    variable = electric_field_x
    sign = negative
    component = x
    block = 'powder_compact graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [electrostatic_calculation_y_graphite_yttria]
    type = PotentialToFieldAux
    gradient_variable = electric_potential
    variable = electric_field_y
    sign = negative
    component = y
    block = 'powder_compact graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [electrostatic_calculation_x_stainless_steel]
    type = PotentialToFieldAux
    gradient_variable = potential_stainless_steel
    variable = electric_field_x
    sign = negative
    component = x
    block = stainless_steel
  []
  [electrostatic_calculation_y_stainless_steel]
    type = PotentialToFieldAux
    gradient_variable = potential_stainless_steel
    variable = electric_field_y
    sign = negative
    component = y
    block = stainless_steel
  []
  [current_density_graphite_yttria]
    type = ADCurrentDensity
    variable = current_density
    potential = electric_potential
    block = 'powder_compact graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [current_density_stainless_steel]
    type = ADCurrentDensity
    variable = current_density
    potential = potential_stainless_steel
    block = stainless_steel
  []
  [microapp_potential]
    type = ParsedAux
    variable = microapp_potential
    coupled_variables = electric_potential
    expression = 'electric_potential*1e9' #convert from V to nV
    block = 'powder_compact'
  []
  [E_x]
    type = VariableGradientComponent
    variable = E_x
    gradient_variable = electric_potential
    component = x
    block = 'powder_compact'
  []
  [E_y]
    type = VariableGradientComponent
    variable = E_y
    gradient_variable = electric_potential
    component = y
    block = 'powder_compact'
  []
  [yttria_current_density_forBC_microapp]
    type = ParsedAux
    variable = yttria_current_density_forBC_microapp
    coupled_variables = 'electrical_conductivity E_y'
    expression = '-1.0*electrical_conductivity*E_y'
    block = 'powder_compact'
  []

  [heat_flux_graphite_x]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature
    variable = heatflux_graphite_x
    component = x
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [heat_flux_graphite_y]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature
    variable = heatflux_graphite_y
    component = y
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [heat_flux_stainless_x]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature_stainless_steel
    variable = heatflux_stainless_steel_x
    block = stainless_steel
    component = x
  []
  [heat_flux_stainless_y]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature_stainless_steel
    variable = heatflux_stainless_steel_y
    block = stainless_steel
    component = y
  []
  [heat_flux_yttria_x]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature
    variable = heatflux_yttria_x
    component = x
    block = 'powder_compact'
  []
  [heat_flux_yttria_y]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature
    variable = heatflux_yttria_y
    component = y
    block = 'powder_compact'
  []

  [heat_transfer_radiation_graphite]
    type = ParsedAux
    variable = heat_transfer_radiation
    boundary = 'outer_radiative_spacers outer_die_wall radiative_upper_plunger radiative_lower_plunger'
    coupled_variables = 'temperature'
    constant_names = 'boltzmann epsilon temperature_farfield' #published emissivity for graphite is 0.85
    constant_expressions = '5.67e-8 0.85 293.0' #roughly room temperature, which is probably too cold
    expression = '-boltzmann*epsilon*(temperature^4-temperature_farfield^4)'
  []
  [heat_transfer_radiation_stainless_steel]
    type = ParsedAux
    variable = heat_transfer_radiation
    boundary = 'outer_radiative_stainless_steel'
    coupled_variables = 'temperature_stainless_steel'
    constant_names = 'boltzmann epsilon temperature_farfield' #published emissivity for graphite is 0.85
    constant_expressions = '5.67e-8 0.4 293.0' #roughly room temperature, which is probably too cold
    expression = '-boltzmann*epsilon*(temperature_stainless_steel^4-temperature_farfield^4)'
  []
  [thermal_conductivity_graphite]
    type = ADMaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [thermal_conductivity_steel]
    type = ADMaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = 'stainless_steel'
  []
  [thermal_conductivity_yttria]
    type = ADMaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = 'powder_compact'
  []
  [electrical_conductivity_graphite]
    type = ADMaterialRealAux
    variable = electrical_conductivity
    property = electrical_conductivity
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [electrical_conductivity_steel]
    type = ADMaterialRealAux
    variable = electrical_conductivity
    property = electrical_conductivity
    block = 'stainless_steel'
  []
  [electrical_conductivity_yttria]
    type = ADMaterialRealAux
    variable = electrical_conductivity
    property = electrical_conductivity
    block = 'powder_compact'
  []
[]

[BCs]
  [external_surface_stainless]
    type = CoupledVarNeumannBC
    boundary = 'outer_radiative_stainless_steel'
    variable = temperature_stainless_steel
    v = heat_transfer_radiation
  []
  [external_surface_graphite]
    type = CoupledVarNeumannBC
    boundary = 'outer_radiative_spacers outer_die_wall radiative_upper_plunger radiative_lower_plunger'
    variable = temperature
    v = heat_transfer_radiation
  []
  [water_channel]
    type = CoupledConvectiveHeatFluxBC
    boundary = water_channel
    variable = temperature_stainless_steel
    T_infinity = 300 #units of K
    htc = 4725
  []
  [temperature_ram_extremes]
    type = ADDirichletBC
    variable = temperature_stainless_steel
    boundary = 'top_upper_steel_ram bottom_lower_steel_ram'
    value = 293
  []
  [electric_potential_top]
    type = ADFunctionNeumannBC
    variable = potential_stainless_steel
    boundary = top_upper_steel_ram
    ## This will need to be updated to match the Dr. Sinter geometry
    function = 'if(t < 31, (980 / (pi * 0.00155))*((0.141301/4.3625)*t), if(t > 600, 0, 980 / (pi * 0.00155)))' # RMS Current / Cross-sectional Area. Ramping for t < 31s approximately reflects Cincotti et al (DOI: 10.1002/aic.11102) Figure 21(b)
  []
  [electric_potential_bottom]
    type = ADDirichletBC
    variable = potential_stainless_steel
    boundary = bottom_lower_steel_ram
    value = 0
  []
[]

[Functions]
  [graphite_thermal_conductivity_fcn]
    type = ParsedFunction
    expression = '-4.418e-12*t^4+2.904e-8*t^3-4.688e-5*t^2-0.0316*t+119.659'
  []
  # [yttria_thermal_conductivity_fcn] #from the multiapp
  #   type = ParsedFunction
  #   expression = '3214.46/(t-147.73)'
  # []
  [harmonic_mean_thermal_conductivity]
    type = ParsedFunction
    expression = '2*(-4.418e-12*t^4+2.904e-8*t^3-4.688e-5*t^2-0.0316*t+119.659)*(3214.46/(t-147.73))/((-4.418e-12*t^4+2.904e-8*t^3-4.688e-5*t^2-0.0316*t+119.659)+(3214.46/(t-147.73)))'
    # symbol_names = 'k_graphite k_yttria'
    # symbol_values = 'graphite_thermal_conductivity_fcn yttria_thermal_conductivity_fcn'
  []

  [graphite_electrical_conductivity_fcn]
    type = ParsedFunction
    expression = '1.0/(-2.705e-15*t^3+1.263e-11*t^2-1.836e-8*t+1.813e-5)'
  []
  # [electrical_conductivity_fcn]
  #   type = ParsedFunction
  #   # symbol_names = porosity
  #   # symbol_values = initial_porosity
  #   expression = '(1-0.62)*2.0025e4*exp(-1.61/8.617343e-5/t)'
  # []
  [harmonic_mean_electrical_conductivity]
    type = ParsedFunction
    expression = '2*(1.0/(-2.705e-15*t^3+1.263e-11*t^2-1.836e-8*t+1.813e-5))*((1)*2.0025e4*exp(-1.61/8.617343e-5/t))/((1.0/(-2.705e-15*t^3+1.263e-11*t^2-1.836e-8*t+1.813e-5))+((1)*2.0025e4*exp(-1.61/8.617343e-5/t)))'
    # symbol_names = 'k_graphite k_yttria'
    # symbol_values = 'graphite_thermal_conductivity_fcn yttria_thermal_conductivity_fcn'
  []
  [mechanical_pressure_func]
    type = ParsedFunction
    symbol_names = 'radius coolant_radius force'
    symbol_values = '0.04 7.071e-3 6.5' # 'm kN'
    expression = 'force * 1e3 / (pi * (radius^2 - coolant_radius^2))' # (N / m^2)
  []
[]

[InterfaceKernels]
  [electric_contact_conductance_ssg]
    type = ElectrostaticContactCondition
    variable = potential_stainless_steel
    neighbor_var = electric_potential
    boundary = 'upper_ram_spacer_interface lower_ram_spacer_interface'
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
  []
  [thermal_contact_conductance_calculated_ssg]
    type = ThermalContactCondition
    variable = temperature_stainless_steel
    neighbor_var = temperature
    primary_potential = potential_stainless_steel
    secondary_potential = electric_potential
    primary_thermal_conductivity = thermal_conductivity
    secondary_thermal_conductivity = thermal_conductivity
    primary_electrical_conductivity = electrical_conductivity
    secondary_electrical_conductivity = electrical_conductivity
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
    boundary = 'upper_ram_spacer_interface lower_ram_spacer_interface'
  []
  [electric_contact_conductance_gss]
    type = ElectrostaticContactCondition
    variable = electric_potential
    neighbor_var = potential_stainless_steel
    boundary = 'upper_ram_spacer_interface lower_ram_spacer_interface'
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
  []
  [thermal_contact_conductance_calculated_gss]
    type = ThermalContactCondition
    variable = temperature
    neighbor_var = temperature_stainless_steel
    primary_potential = electric_potential
    secondary_potential = potential_stainless_steel
    primary_thermal_conductivity = thermal_conductivity
    secondary_thermal_conductivity = thermal_conductivity
    primary_electrical_conductivity = electrical_conductivity
    secondary_electrical_conductivity = electrical_conductivity
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
    boundary = 'upper_ram_spacer_interface lower_ram_spacer_interface'
  []
[]

[ThermalContact]
  [upper_plunger_powder_electric]
    type = GapHeatTransfer
    primary = bottom_upper_plunger
    secondary = top_powder_compact
    variable = electric_potential
    quadrature = true
    emissivity_primary = 0.0 #not applicable for electric potential
    emissivity_secondary = 0.0 #not applicable for electric potential
    gap_geometry_type = PLATE
    # gap_conductivity = 1e-23 #rough harmonic mean
    gap_conductivity_function = 'harmonic_mean_electrical_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [powder_bottom_plunger_electric]
    type = GapHeatTransfer
    primary = top_lower_plunger
    secondary = bottom_powder_compact
    variable = electric_potential
    quadrature = true
    emissivity_primary = 0.0 #not applicable for electric potential
    emissivity_secondary = 0.0 #not applicable for electric potential
    gap_geometry_type = PLATE
    # gap_conductivity = 1e-23 #rough harmonic mean
    gap_conductivity_function = 'harmonic_mean_electrical_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [powder_die_electric]
    type = GapHeatTransfer
    primary = inner_die_wall
    secondary = outer_powder_compact
    variable = electric_potential
    quadrature = true
    emissivity_primary = 0.0 #0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.0 #0.85
    # gap_geometry_type = PLATE # Not for vertical surfaces
    # gap_conductivity = 1e-23 #rough harmonic mean
    gap_conductivity_function = 'harmonic_mean_electrical_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [combined_plungers_die_electric]
    type = GapHeatTransfer
    primary = inner_die_wall
    secondary = combined_outer_plungers
    variable = electric_potential
    quadrature = true
    emissivity_primary = 0.0 #not applicable for electric potential
    emissivity_secondary = 0.0 #not applicable for electric potential
    # gap_geometry_type = PLATE # Not for vertical surfaces
    # gap_conductivity = 8.5e4 #graphite, at 500K
    gap_conductivity_function = 'graphite_electrical_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []

##Thermal Contact between gapped graphite die components
  [upper_plunger_spacer_gap_thermal]
    type = GapHeatTransfer
    primary = spacer_facing_upper_plunger
    secondary = plunger_facing_upper_spacer
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    gap_geometry_type = PLATE
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [upper_diewall_spacer_gap_thermal]
    type = GapHeatTransfer
    primary = top_die_wall
    secondary = bottom_upper_spacer
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    gap_geometry_type = PLATE
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [lower_plunger_spacer_gap_thermal]
    type = GapHeatTransfer
    primary = spacer_facing_lower_plunger
    secondary = plunger_facing_lower_spacer
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    gap_geometry_type = PLATE
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [lower_diewall_spacer_gap_thermal]
    type = GapHeatTransfer
    primary = bottom_die_wall
    secondary = top_lower_spacer
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    gap_geometry_type = PLATE
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [upper_plunger_diewall_gap_thermal]
    type = GapHeatTransfer
    primary = inner_die_wall ### paired temperature doesn't show on inner die wall, but temperature profile looks reasonable
    secondary = die_wall_facing_upper_plunger
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    # gap_geometry_type = PLATE # Not for vertical surfaces
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [lower_plunger_diewall_gap_thermal]
    type = GapHeatTransfer
    primary = inner_die_wall ### paired temperature doesn't show on inner die wall, but temperature profile looks reasonable
    secondary = die_wall_facing_lower_plunger
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    # gap_geometry_type = PLATE # Not for vertical surfaces
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []

## Thermal Contact between touching components of powder and die
  [upper_plunger_powder_thermal]
    type = GapHeatTransfer
    primary = bottom_upper_plunger
    secondary = top_powder_compact
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.3 #estimated from McMahon and Wilder, High Temperature Spectral Emissivity of Yttrium, Samarium, Gadolinium, Ebrium and Lutetium Oxides (1963) Atomic Energy Commission, IS-578, Figure 12
    gap_geometry_type = PLATE
    gap_conductivity_function = 'harmonic_mean_thermal_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [powder_bottom_plunger_thermal]
    type = GapHeatTransfer
    primary = top_lower_plunger
    secondary = bottom_powder_compact #expect more heat transfer from the die to the powder
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.3 #estimated from McMahon and Wilder, High Temperature Spectral Emissivity of Yttrium, Samarium, Gadolinium, Ebrium and Lutetium Oxides (1963) Atomic Energy Commission, IS-578, Figure 12
    gap_geometry_type = PLATE
    gap_conductivity_function = 'harmonic_mean_thermal_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [powder_die_thermal]
    type = GapHeatTransfer
    primary = inner_die_wall
    secondary = outer_powder_compact
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.3 #estimated from McMahon and Wilder, High Temperature Spectral Emissivity of Yttrium, Samarium, Gadolinium, Ebrium and Lutetium Oxides (1963) Atomic Energy Commission, IS-578, Figure 12
    # gap_geometry_type = PLATE # Not for vertical surfaces
    gap_conductivity_function = 'harmonic_mean_thermal_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [combined_plungers_die_thermal]
    type = GapHeatTransfer
    primary = inner_die_wall
    secondary = combined_outer_plungers
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    # gap_geometry_type = PLATE # Not for vertical surfaces
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
[]

[Materials]
  ## graphite blocks
  [graphite_thermal]
    type = ADGraphiteThermal
    temperature = temperature
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [graphite_density]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = 1.750e3 #in kg/m^3 from Cincotti et al 2007, Table 2, doi:10.1002/aic
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [graphite_electrical_conductivity]
    type = ADGraphiteElectricalConductivity
    temperature = temperature
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []

  #stainless_steel
  [stainless_steel_thermal]
    type = ADStainlessSteelThermal
    temperature = temperature_stainless_steel
    block = stainless_steel
  []
  [stainless_steel_density]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = 8000
    block = stainless_steel
  []
  [stainless_steel_electrical_conductivity]
    type = ADStainlessSteelElectricalConductivity
    temperature = temperature_stainless_steel
    block = stainless_steel
  []

  # harmonic mean of graphite and stainless steel hardness
  [mean_hardness_graphite_stainless_steel]
    type = GraphiteStainlessMeanHardness
    block = 'graphite_spacers upper_plunger lower_plunger die_wall stainless_steel'
  []

  ## yttria powder compact
  [yttria_thermal_conductivity]
    type = ADParsedMaterial
    # coupled_variables = 'temperature'
    # expression = '3214.46 / (temperature - 147.73)' #in W/(m-K) #Given from Larry's curve fitting, data from Klein and Croft, JAP, v. 38, p. 1603 and UC report "For Computer Heat Conduction Calculations - A compilation of thermal properties data" by A.L. Edwards, UCRL-50589 (1969)
    coupled_variables = 'yttria_thermal_conductivity_aeh'
    expression = 'yttria_thermal_conductivity_aeh' #in W/(m-K) directly, for now
    property_name = 'thermal_conductivity'
    output_properties = thermal_conductivity
    outputs = 'exodus'
    block = powder_compact
  []
  [yttria_specific_heat_capacity]
    type = ADParsedMaterial
    property_name = heat_capacity
    coupled_variables = 'yttria_heat_capacity_volume_avg'
    expression = 'yttria_heat_capacity_volume_avg' #in J/(K-kg)
    # output_properties = yttria_specific_heat_capacity
    # outputs = 'exodus'
    block = powder_compact
  []
  [yttria_density]
    type = ADParsedMaterial
    property_name = 'density'
    coupled_variables = 'yttria_density_volume_avg'
    expression = 'yttria_density_volume_avg'
    # output_properties = yttria_density
    # outputs = 'exodus'
    block = powder_compact
  []
  [electrical_conductivity]
    type = ADParsedMaterial
    coupled_variables = 'yttria_sigma_aeh'
    expression = 'yttria_sigma_aeh*1.602e-10' #converts from eV/(V^2 s nm) to units of J/(V^2-m-s)
    property_name = 'electrical_conductivity'
    output_properties = electrical_conductivity
    outputs = 'exodus'
    block = powder_compact
    # type = ADDerivativeParsedMaterial
    # property_name = electrical_conductivity
    # coupled_variables = 'temperature'
    # constant_names =       'Q_elec  kB            prefactor_solid  initial_porosity'
    # constant_expressions = '1.61    8.617343e-5        1.25e-4           0.38'
    # expression = '(1-initial_porosity) * prefactor_solid * exp(-Q_elec/kB/temperature) * 1.602e8' # in eV/(nV^2 s nm) per chat with Larry, last term converts to units of J/(V^2-m-s)
  []

  # Material property converter for DiffusionFluxAux object
  [converter]
    type = MaterialADConverter
    ad_props_in = thermal_conductivity
    reg_props_out = nonad_thermal_conductivity
  []

  [Q_SI]
    type = ParsedMaterial
    property_name = Q_SI
    coupled_variables = 'Q_from_sub'
    expression = 'Q_from_sub / 260 / 260 * 0.1602' #divide by domain size and unit conversion to go from eV/S to J/m^3/s
    outputs = 'exodus'
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  automatic_scaling = true
  line_search = 'none'
  # compute_scaling_once = false

  # force running options
  # petsc_options_iname = '-pc_type -snes_linesearch_type -pc_factor_shift_type -pc_factor_shift_amount'
  # petsc_options_value = 'lu       basic                 NONZERO               1e-15'

  # #mechanical contact options
  # petsc_options = '-snes_ksp_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = ' lu       superlu_dist'
  petsc_options = '-snes_converged_reason -ksp_converged_reason'

  nl_forced_its = 1
  nl_rel_tol = 2e-5 #1e-6 #2e-5 for with mechanics #was 1e-10, for temperature only
  nl_abs_tol = 2e-12 #was 1e-12
  nl_max_its = 20
  l_max_its = 50
  dtmin = 1.0e-4

  end_time = 900 #600 #900 #15 minutes, rule of thumb from Dennis is 10 minutes
  [Quadrature]
    order = FIFTH #required for thermal and mechanical node-face contact
    side_order = SEVENTH
  []
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 0.05
    optimal_iterations = 8
    iteration_window = 2
  []
[]

[Postprocessors]
  [temperature_pp]
    type = AverageNodalVariableValue
    variable = temperature
    block = 'graphite_spacers upper_plunger lower_plunger die_wall powder_compact'
  []
  [temperature_stainless_steel_pp]
    type = AverageNodalVariableValue
    variable = temperature_stainless_steel
    block = stainless_steel
  []
  [graphite_thermal_conductivity]
    type = ElementAverageValue
    variable = thermal_conductivity
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [graphite_electrical_conductivity]
    type = ElementAverageValue
    variable = electrical_conductivity
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [steel_thermal_conductivity]
    type = ElementAverageValue
    variable = thermal_conductivity
    block = stainless_steel
  []
  [steel_electrical_conductivity]
    type = ElementAverageValue
    variable = electrical_conductivity
    block = stainless_steel
  []
  [yttria_thermal_conductivity]
    type = ElementAverageValue
    variable = thermal_conductivity
    block = powder_compact
  []
  [yttria_electrical_conductivity]
    type = ElementAverageValue
    variable = electrical_conductivity
    block = powder_compact
  []
[]

[MultiApps]
  [micro]
    type = TransientMultiApp
    # type = CentroidMultiApp # lauches one in the middle of each element so don't need to give positions
    #can specify the number of procs
    max_procs_per_app = 10 #Can specify more than this, but it will degrade convergence
    app_type = MalamuteApp
    positions = '0.00125 0.034 0' #roughly the center of element 117 in this mesh
    input_files = micro_yttria_thermoelectric_twoway_lots_controls.i
    catch_up = true
    execute_on = TIMESTEP_BEGIN #the default
  []
[]

[Transfers]
  [keff_from_sub]
    type = MultiAppPostprocessorInterpolationTransfer
    from_multi_app = micro
    variable = yttria_thermal_conductivity_aeh
    power = 2 #2 is the default value, tutorial uses 1
    postprocessor = k_AEH_average
  []
  [sigma_aeh_eff_from_sub]
    type = MultiAppPostprocessorInterpolationTransfer
    from_multi_app = micro
    variable = yttria_sigma_aeh
    power = 2 #2 is the default value, tutorial uses 1
    postprocessor = sigma_y_AEH
  []
  [Q_from_sub]
    type = MultiAppPostprocessorInterpolationTransfer
    from_multi_app = micro
    variable = Q_from_sub #This is the integrated heat produced in the phase-field simulation in eV/m/s
    power = 2 #2 is the default value, tutorial uses 1
    postprocessor = Q_joule_total
  []
  [temperature_to_sub]
    type = MultiAppVariableValueSampleTransfer
    to_multi_app = micro
    source_variable = temperature
    variable = T
  []
  [temperature_to_sub_postproc]
    type = MultiAppVariableValueSamplePostprocessorTransfer
    to_multi_app = micro
    source_variable = temperature
    postprocessor = T_postproc
  []
  [potential_to_sub_postproc]
    type = MultiAppVariableValueSamplePostprocessorTransfer
    to_multi_app = micro
    source_variable = electric_potential
    postprocessor = V_postproc
  []
  [micro_field_pp_to_sub]
    type = MultiAppVariableValueSamplePostprocessorTransfer
    to_multi_app = micro
    source_variable = E_y
    postprocessor = Ey_in
  []
[]

[Outputs]
  csv = true
  exodus = true
  perf_graph = true
  # [ckpt]
  #   type =Checkpoint
  #   time_step_interval = 1
  #   num_files = 2
  # []
[]

(examples/sps/multiapp/complete_geometry/electrothermal/engineering_scale_electrothermal_goofycontact.i)

# Approximate Dr. Sinter Geometry
# Configuration: Two stainless steel electrodes (with cooling channels) with two
#                large spacers, and one die-and-plungers assembly.
#                RZ, axisymmetric.
# BCs:
#    Potential:
#       V (top electrode, top surface) --> Neumann condition specifiying potential
#                                          based on applied RMS Current (980 A)
#                                          and cross-sectional electrode area (see
#                                          elec_top BC below). Current turned off
#                                          at 1200 s.
#       V (bottom electrode, bottom surface) = 0 V
#       V (elsewhere) --> natural boundary conditions (no current external to circuit)
#    Temperature:
#       T (top electrode, top surface) = 293 K
#       T (bottom electrode, bottom surface) = 293 K
#       T (vertical right side) --> simple radiative BC into black body at 293 K
#       T (horizonal surfaces) --> GapHeatTransfer
#       T (water channel) --> simple convective BC into "water" with heat transfer coefficient parameter from Cincotti and temperature of 293 K
#       T (elsewhere) --> natural boundary conditions
#    Mechanics (not added yet):
#       Displacement: --> pinned axial centerpoint along interior surface of powder block, outer die wall
#       Pressure: top and bottom surfaces of upper and lower rams, 20.7 MPa
#                 ramped quickly at simulation start, then held constant
# Interface Conditions:
#       V (SS-G upper and G-SS lower interface) --> ElectrostaticContactCondition
#       T (SS-G upper and G-SS lower interface) --> ThermalContactCondition
#       V (graphite-yttria and graphite-graphite) --> GapHeatTransfer
#       T (graphite-yttria and graphite-graphite) --> GapHeatTransfer
#       P (graphite-yttria and graphite-graphite) --> Mechanical Contact
# Initial Conditions:
#       V = default (0 V)
#       T = 293 K
# Applied Mechanical Load: needs to be updated to 6.5kN, was 3.5 kN
#
# Reference for graphite, stainless steel: Cincotti et al, DOI 10.1002/aic.11102
# Assorted references for yttria, listed as comments in input file

initial_temperature = 293 #roughly 600C where the pyrometer kicks in
#initial_porosity=0.36 #Maximum random jammed packing, Donev et al (2004) Science Magazine

[Mesh]
  file = ../../../meshes/drsinter_nodalcontact_2d.e
  coord_type = RZ
  construct_side_list_from_node_list = true
  patch_update_strategy = iteration
  patch_size = 20
  # ghosting_patch_size = 5*patch_size
[]

[Problem]
  type = ReferenceResidualProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
[]

[Variables]
  [temperature_stainless_steel]
    initial_condition = ${initial_temperature}
    block = stainless_steel
  []
  [temperature]
    initial_condition = ${initial_temperature}
    block = 'graphite_spacers lower_plunger upper_plunger die_wall powder_compact'
  []
  [potential_stainless_steel]
    block = stainless_steel
  []
  [electric_potential]
    block = 'graphite_spacers lower_plunger upper_plunger die_wall powder_compact'
  []
[]

[AuxVariables]
  [electric_field_x]
    family = MONOMIAL #prettier pictures with smoother values
    order = FIRST
  []
  [electric_field_y]
    family = MONOMIAL
    order = FIRST
  []
  [current_density]
    family = NEDELEC_ONE
    order = FIRST
  []
  [yttria_sigma_aeh]
    initial_condition = 2.0e-10 #in units eV/((nV)^2-s-nm)
    block = 'powder_compact'
  []
  [microapp_potential]
    #converted to microapp electronVolts units
    block = 'powder_compact'
  []
  [E_x]
    order = FIRST
    family = MONOMIAL
    block = 'powder_compact'
  []
  [E_y]
    order = FIRST
    family = MONOMIAL
    block = 'powder_compact'
  []
  [yttria_current_density_forBC_microapp]
    order = FIRST
    family = MONOMIAL
    block = 'powder_compact'
  []

  # [T_infinity]
  #   initial_condition = ${initial_temperature}
  # []
  [heatflux_graphite_x]
    family = MONOMIAL
    order = CONSTANT
    block = 'graphite_spacers lower_plunger upper_plunger die_wall'
  []
  [heatflux_graphite_y]
    family = MONOMIAL
    order = CONSTANT
    block = 'graphite_spacers lower_plunger upper_plunger die_wall'
  []
  [heatflux_stainless_steel_x]
    family = MONOMIAL
    order = CONSTANT
    block = stainless_steel
  []
  [heatflux_stainless_steel_y]
    family = MONOMIAL
    order = CONSTANT
    block = stainless_steel
  []
  [heatflux_yttria_x]
    family = MONOMIAL
    order = CONSTANT
    block = 'powder_compact'
  []
  [heatflux_yttria_y]
    family = MONOMIAL
    order = CONSTANT
    block = 'powder_compact'
  []
  [yttria_thermal_conductivity_aeh]
    initial_condition = 0.4
  []
  [yttria_heat_capacity_volume_avg]
    initial_condition = 447.281 #842.2 at 1600K probably? # at 1500K #568.73 at 1000K #447.281 # at 293K
  []
  [yttria_density_volume_avg]
    initial_condition = 3206.4 ##5010.0*(1-${initial_porosity}) #in kg/m^3
  []
  [heat_transfer_radiation]
  []
  # [thermal_conductivity]
  #   family = MONOMIAL
  #   order = FIRST
  # []
[]

[Kernels]
  [HeatDiff_graphite]
    type = ADHeatConduction
    variable = temperature
    thermal_conductivity = thermal_conductivity
    extra_vector_tags = 'ref'
    block = 'graphite_spacers lower_plunger upper_plunger die_wall'
  []
  [HeatTdot_graphite]
    type = ADHeatConductionTimeDerivative
    variable = temperature
    specific_heat = heat_capacity
    density_name = density
    extra_vector_tags = 'ref'
    block = 'graphite_spacers lower_plunger upper_plunger die_wall'
  []
  [JouleHeating_graphite]
    type = ADJouleHeatingSource
    variable = temperature
    elec = electric_potential
    electrical_conductivity = electrical_conductivity
    use_displaced_mesh = true
    extra_vector_tags = 'ref'
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [electric_graphite]
    type = ADMatDiffusion
    variable = electric_potential
    diffusivity = electrical_conductivity
    extra_vector_tags = 'ref'
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []

  [HeatDiff_stainless_steel]
    type = ADHeatConduction
    variable = temperature_stainless_steel
    thermal_conductivity = thermal_conductivity
    extra_vector_tags = 'ref'
    block = stainless_steel
  []
  [HeatTdot_stainless_steel]
    type = ADHeatConductionTimeDerivative
    variable = temperature_stainless_steel
    specific_heat = heat_capacity
    density_name = density
    extra_vector_tags = 'ref'
    block = stainless_steel
  []
  [HeatSource_stainless_steel]
    type = ADJouleHeatingSource
    variable = temperature_stainless_steel
    elec = potential_stainless_steel
    electrical_conductivity = electrical_conductivity
    use_displaced_mesh = true
    extra_vector_tags = 'ref'
    block = stainless_steel
  []
  [electric_stainless_steel]
    type = ADMatDiffusion
    variable = potential_stainless_steel
    diffusivity = electrical_conductivity
    extra_vector_tags = 'ref'
    block = stainless_steel
  []

  [HeatDiff_yttria]
    type = ADHeatConduction
    variable = temperature
    thermal_conductivity = thermal_conductivity #use parsed material property, hope it works
    extra_vector_tags = 'ref'
    block = powder_compact
  []
  [HeatTdot_yttria]
    type = ADHeatConductionTimeDerivative
    variable = temperature
    specific_heat = heat_capacity #use parsed material property
    density_name = density
    extra_vector_tags = 'ref'
    block = powder_compact
  []
  [JouleHeating_yttria]
    type = ADJouleHeatingSource
    variable = temperature
    elec = electric_potential
    electrical_conductivity = electrical_conductivity
    use_displaced_mesh = true
    extra_vector_tags = 'ref'
    block = powder_compact
  []
  [electric_yttria]
    type = ADMatDiffusion
    variable = electric_potential
    diffusivity = electrical_conductivity
    extra_vector_tags = 'ref'
    block = powder_compact
  []
[]

[AuxKernels]
  [electrostatic_calculation_x_graphite_yttria]
    type = PotentialToFieldAux
    gradient_variable = electric_potential
    variable = electric_field_x
    sign = negative
    component = x
    block = 'powder_compact graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [electrostatic_calculation_y_graphite_yttria]
    type = PotentialToFieldAux
    gradient_variable = electric_potential
    variable = electric_field_y
    sign = negative
    component = y
    block = 'powder_compact graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [electrostatic_calculation_x_stainless_steel]
    type = PotentialToFieldAux
    gradient_variable = potential_stainless_steel
    variable = electric_field_x
    sign = negative
    component = x
    block = stainless_steel
  []
  [electrostatic_calculation_y_stainless_steel]
    type = PotentialToFieldAux
    gradient_variable = potential_stainless_steel
    variable = electric_field_y
    sign = negative
    component = y
    block = stainless_steel
  []
  [current_density_graphite_yttria]
    type = ADCurrentDensity
    variable = current_density
    potential = electric_potential
    block = 'powder_compact graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [current_density_stainless_steel]
    type = ADCurrentDensity
    variable = current_density
    potential = potential_stainless_steel
    block = stainless_steel
  []
  [microapp_potential]
    type = ParsedAux
    variable = microapp_potential
    coupled_variables = electric_potential
    expression = 'electric_potential*1e9' #convert from V to nV
    block = 'powder_compact'
  []
  [E_x]
    type = VariableGradientComponent
    variable = E_x
    gradient_variable = electric_potential
    component = x
    block = 'powder_compact'
  []
  [E_y]
    type = VariableGradientComponent
    variable = E_y
    gradient_variable = electric_potential
    component = y
    block = 'powder_compact'
  []
  [yttria_current_density_forBC_microapp]
    type = ParsedAux
    variable = yttria_current_density_forBC_microapp
    coupled_variables = 'electrical_conductivity E_y'
    expression = '-1.0*electrical_conductivity*E_y'
    block = 'powder_compact'
  []

  [heat_flux_graphite_x]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature
    variable = heatflux_graphite_x
    component = x
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [heat_flux_graphite_y]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature
    variable = heatflux_graphite_y
    component = y
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [heat_flux_stainless_x]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature_stainless_steel
    variable = heatflux_stainless_steel_x
    block = stainless_steel
    component = x
  []
  [heat_flux_stainless_y]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature_stainless_steel
    variable = heatflux_stainless_steel_y
    block = stainless_steel
    component = y
  []
  [heat_flux_yttria_x]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature
    variable = heatflux_yttria_x
    component = x
    block = 'powder_compact'
  []
  [heat_flux_yttria_y]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature
    variable = heatflux_yttria_y
    component = y
    block = 'powder_compact'
  []

  [heat_transfer_radiation_graphite]
    type = ParsedAux
    variable = heat_transfer_radiation
    boundary = 'outer_radiative_spacers outer_die_wall radiative_upper_plunger radiative_lower_plunger'
    coupled_variables = 'temperature'
    constant_names = 'boltzmann epsilon temperature_farfield' #published emissivity for graphite is 0.85
    constant_expressions = '5.67e-8 0.85 293.0' #roughly room temperature, which is probably too cold
    expression = '-boltzmann*epsilon*(temperature^4-temperature_farfield^4)'
  []
  [heat_transfer_radiation_stainless_steel]
    type = ParsedAux
    variable = heat_transfer_radiation
    boundary = 'outer_radiative_stainless_steel'
    coupled_variables = 'temperature'
    constant_names = 'boltzmann epsilon temperature_farfield' #published emissivity for graphite is 0.85
    constant_expressions = '5.67e-8 0.4 293.0' #roughly room temperature, which is probably too cold
    expression = '-boltzmann*epsilon*(temperature^4-temperature_farfield^4)'
  []
  [thermal_conductivity_graphite]
    type = ADMaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [thermal_conductivity_steel]
    type = ADMaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = 'stainless_steel'
  []
  [thermal_conductivity_yttria]
    type = ADMaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = 'powder_compact'
  []
  [electrical_conductivity_graphite]
    type = ADMaterialRealAux
    variable = electrical_conductivity
    property = electrical_conductivity
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [electrical_conductivity_steel]
    type = ADMaterialRealAux
    variable = electrical_conductivity
    property = electrical_conductivity
    block = 'stainless_steel'
  []
  [electrical_conductivity_yttria]
    type = ADMaterialRealAux
    variable = electrical_conductivity
    property = electrical_conductivity
    block = 'powder_compact'
  []
[]

[BCs]
  [external_surface_stainless]
    type = CoupledVarNeumannBC
    boundary = 'outer_radiative_stainless_steel'
    variable = temperature_stainless_steel
    v = heat_transfer_radiation
  []
  [external_surface_graphite]
    type = CoupledVarNeumannBC
    boundary = 'outer_radiative_spacers outer_die_wall radiative_upper_plunger radiative_lower_plunger'
    variable = temperature
    v = heat_transfer_radiation
  []
  [water_channel]
    type = CoupledConvectiveHeatFluxBC
    boundary = water_channel
    variable = temperature_stainless_steel
    T_infinity = ${initial_temperature}
    htc = 4725
  []
  [temperature_ram_extremes]
    type = ADDirichletBC
    variable = temperature_stainless_steel
    boundary = 'top_upper_steel_ram bottom_lower_steel_ram'
    value = 293
  []
  [electric_potential_top]
    type = ADFunctionNeumannBC
    variable = potential_stainless_steel
    boundary = top_upper_steel_ram
    ## This will need to be updated to match the Dr. Sinter geometry
    function = 'if(t < 31, (980 / (pi * 0.00155))*((0.141301/4.3625)*t), if(t > 600, 0, 980 / (pi * 0.00155)))' # RMS Current / Cross-sectional Area. Ramping for t < 31s approximately reflects Cincotti et al (DOI: 10.1002/aic.11102) Figure 21(b)
  []
  [electric_potential_bottom]
    type = ADDirichletBC
    variable = potential_stainless_steel
    boundary = bottom_lower_steel_ram
    value = 0
  []
[]

[Functions]
  [graphite_thermal_conductivity_fcn]
    type = ParsedFunction
    expression = '-4.418e-12*t^4+2.904e-8*t^3-4.688e-5*t^2-0.0316*t+119.659'
  []
  # [yttria_thermal_conductivity_fcn] #from the multiapp
  #   type = ParsedFunction
  #   expression = '3214.46/(t-147.73)'
  # []
  [harmonic_mean_thermal_conductivity]
    type = ParsedFunction
    expression = '2*(-4.418e-12*t^4+2.904e-8*t^3-4.688e-5*t^2-0.0316*t+119.659)*(3214.46/(t-147.73))/((-4.418e-12*t^4+2.904e-8*t^3-4.688e-5*t^2-0.0316*t+119.659)+(3214.46/(t-147.73)))'
    # symbol_names = 'k_graphite k_yttria'
    # symbol_values = 'graphite_thermal_conductivity_fcn yttria_thermal_conductivity_fcn'
  []

  [graphite_electrical_conductivity_fcn]
    type = ParsedFunction
    expression = '1.0/(-2.705e-15*t^3+1.263e-11*t^2-1.836e-8*t+1.813e-5)'
  []
  # [electrical_conductivity_fcn]
  #   type = ParsedFunction
  #   # symbol_names = porosity
  #   # symbol_values = initial_porosity
  #   expression = '(1-0.62)*2.0025e4*exp(-1.61/8.617343e-5/t)'
  # []
  [harmonic_mean_electrical_conductivity]
    type = ParsedFunction
    expression = '2*(1.0/(-2.705e-15*t^3+1.263e-11*t^2-1.836e-8*t+1.813e-5))*((1)*2.0025e4*exp(-1.61/8.617343e-5/t))/((1.0/(-2.705e-15*t^3+1.263e-11*t^2-1.836e-8*t+1.813e-5))+((1)*2.0025e4*exp(-1.61/8.617343e-5/t)))'
    # symbol_names = 'k_graphite k_yttria'
    # symbol_values = 'graphite_thermal_conductivity_fcn yttria_thermal_conductivity_fcn'
  []
  [mechanical_pressure_func]
    type = ParsedFunction
    symbol_names = 'radius force'
    symbol_values = '0.04 6.5' # 'm kN'
    expression = 'force * 1e3 / (pi * radius^2)' # (N / m^2)
  []
[]

[InterfaceKernels]
  [electric_contact_conductance_ssg]
    type = ElectrostaticContactCondition
    variable = potential_stainless_steel
    neighbor_var = electric_potential
    boundary = 'upper_ram_spacer_interface lower_ram_spacer_interface'
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
  []
  [thermal_contact_conductance_calculated_ssg]
    type = ThermalContactCondition
    variable = temperature_stainless_steel
    neighbor_var = temperature
    primary_potential = potential_stainless_steel
    secondary_potential = electric_potential
    primary_thermal_conductivity = thermal_conductivity
    secondary_thermal_conductivity = thermal_conductivity
    primary_electrical_conductivity = electrical_conductivity
    secondary_electrical_conductivity = electrical_conductivity
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
    boundary = 'upper_ram_spacer_interface lower_ram_spacer_interface'
  []
  [electric_contact_conductance_gss]
    type = ElectrostaticContactCondition
    variable = electric_potential
    neighbor_var = potential_stainless_steel
    boundary = 'upper_ram_spacer_interface lower_ram_spacer_interface'
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
  []
  [thermal_contact_conductance_calculated_gss]
    type = ThermalContactCondition
    variable = temperature
    neighbor_var = temperature_stainless_steel
    primary_potential = electric_potential
    secondary_potential = potential_stainless_steel
    primary_thermal_conductivity = thermal_conductivity
    secondary_thermal_conductivity = thermal_conductivity
    primary_electrical_conductivity = electrical_conductivity
    secondary_electrical_conductivity = electrical_conductivity
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
    boundary = 'upper_ram_spacer_interface lower_ram_spacer_interface'
  []
[]

[ThermalContact]
  [upper_plunger_powder_electric]
    type = GapHeatTransfer
    primary = bottom_upper_plunger
    secondary = top_powder_compact
    variable = electric_potential
    quadrature = true
    emissivity_primary = 0.0 #not applicable for electric potential
    emissivity_secondary = 0.0 #not applicable for electric potential
    gap_geometry_type = PLATE
    # gap_conductivity = 1e-23 #rough harmonic mean
    gap_conductivity_function = 'harmonic_mean_electrical_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [powder_bottom_plunger_electric]
    type = GapHeatTransfer
    primary = top_lower_plunger
    secondary = bottom_powder_compact
    variable = electric_potential
    quadrature = true
    emissivity_primary = 0.0 #not applicable for electric potential
    emissivity_secondary = 0.0 #not applicable for electric potential
    gap_geometry_type = PLATE
    # gap_conductivity = 1e-23 #rough harmonic mean
    gap_conductivity_function = 'harmonic_mean_electrical_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [combined_plungers_die_electric]
    #to reproduce a more symmetric temperature profile, move this block to the bottom of this ThermalContact section
    type = GapHeatTransfer
    primary = inner_die_wall
    secondary = combined_outer_plungers
    variable = electric_potential
    quadrature = true
    emissivity_primary = 0.0 #not applicable for electric potential
    emissivity_secondary = 0.0 #not applicable for electric potential
    # gap_geometry_type = PLATE # Not for vertical surfaces
    # gap_conductivity = 8.5e4 #graphite, at 500K
    gap_conductivity_function = 'graphite_electrical_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
    # appended_property_name = plungers
  []
  [powder_die_electric]
    type = GapHeatTransfer
    primary = inner_die_wall
    secondary = outer_powder_compact
    variable = electric_potential
    quadrature = true
    emissivity_primary = 0.0 #0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.0 #0.85
    # gap_geometry_type = PLATE # Not for vertical surfaces
    # gap_conductivity = 1e-23 #rough harmonic mean
    gap_conductivity_function = 'harmonic_mean_electrical_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
    # appended_property_name = powder
  []
[]

# ## Thermal Contact between touching components of powder and die
# [ThermalContact]
# [lower_plunger_die_thermal]
#   type = GapHeatTransfer
#   primary = inner_die_wall
#   secondary = outer_lower_plunger
#   variable = temperature
#   quadrature = true
#   emissivity_primary = 0.85 #cincotti 2007, table 2
#   emissivity_secondary = 0.85
#   # gap_geometry_type = PLATE # Not for vertical surfaces
#   gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
#   gap_conductivity_function_variable = temperature
#   normal_smoothing_distance = 0.1
# []
# [top_plunger_die_thermal]
#   type = GapHeatTransfer
#   primary = inner_die_wall
#   secondary = outer_upper_plunger
#   variable = temperature
#   quadrature = true
#   emissivity_primary = 0.85 #cincotti 2007, table 2
#   emissivity_secondary = 0.85
#   # gap_geometry_type = PLATE # Not for vertical surfaces
#   gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
#   gap_conductivity_function_variable = temperature
#   normal_smoothing_distance = 0.1
# []
# [upper_plunger_powder_thermal]
#   type = GapHeatTransfer
#   primary = bottom_upper_plunger
#   secondary = top_powder_compact
#   variable = temperature
#   quadrature = true
#   emissivity_primary = 0.85 #cincotti 2007, table 2
#   emissivity_secondary = 0.85
#   gap_geometry_type = PLATE
#   gap_conductivity_function = 'harmonic_mean_thermal_conductivity'
#   gap_conductivity_function_variable = temperature
#   normal_smoothing_distance = 0.1
# []
# [powder_die_thermal]
#   type = GapHeatTransfer
#   primary = inner_die_wall
#   secondary = outer_powder_compact
#   variable = temperature
#   quadrature = true
#   emissivity_primary = 0.85 #cincotti 2007, table 2
#   emissivity_secondary = 0.3 #estimated from McMahon and Wilder, High Temperature Spectral Emissivity of Yttrium, Samarium, Gadolinium, Ebrium and Lutetium Oxides (1963) Atomic Energy Commission, IS-578, Figure 12
#   # gap_geometry_type = PLATE # Not for vertical surfaces
#   gap_conductivity_function = 'harmonic_mean_thermal_conductivity'
#   gap_conductivity_function_variable = temperature
#   normal_smoothing_distance = 0.1
# []
# [powder_bottom_plunger_thermal]
#   type = GapHeatTransfer
#   primary = top_lower_plunger
#   secondary = bottom_powder_compact #expect more heat transfer from the die to the powder
#   variable = temperature
#   quadrature = true
#   emissivity_primary = 0.85 #cincotti 2007, table 2
#   emissivity_secondary = 0.3 #estimated from McMahon and Wilder, High Temperature Spectral Emissivity of Yttrium, Samarium, Gadolinium, Ebrium and Lutetium Oxides (1963) Atomic Energy Commission, IS-578, Figure 12
#   gap_geometry_type = PLATE
#   gap_conductivity_function = 'harmonic_mean_thermal_conductivity'
#   gap_conductivity_function_variable = temperature
#   normal_smoothing_distance = 0.1
# []

# ##Thermal Contact between gapped graphite die components
# [lower_plunger_spacer_gap_thermal]
#   type = GapHeatTransfer
#   primary = spacer_facing_lower_plunger
#   secondary = plunger_facing_lower_spacer
#   variable = temperature
#   quadrature = true
#   emissivity_primary = 0.85 #cincotti 2007, table 2
#   emissivity_secondary = 0.85
#   gap_geometry_type = PLATE
#   gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
#   gap_conductivity_function_variable = temperature
#   normal_smoothing_distance = 0.1
# []
# [lower_plunger_diewall_gap_thermal]
#   type = GapHeatTransfer
#   primary = die_wall_facing_lower_plunger
#   secondary = inner_die_wall
#   variable = temperature
#   quadrature = true
#   emissivity_primary = 0.85 #cincotti 2007, table 2
#   emissivity_secondary = 0.85
#   # gap_geometry_type = PLATE # Not for vertical surfaces
#   gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
#   gap_conductivity_function_variable = temperature
#   normal_smoothing_distance = 0.1
# []
# [lower_diewall_spacer_gap_thermal]
#   type = GapHeatTransfer
#   primary = bottom_die_wall
#   secondary = top_lower_spacer
#   variable = temperature
#   quadrature = true
#   emissivity_primary = 0.85 #cincotti 2007, table 2
#   emissivity_secondary = 0.85
#   gap_geometry_type = PLATE
#   gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
#   gap_conductivity_function_variable = temperature
#   normal_smoothing_distance = 0.1
# []
# [upper_plunger_spacer_gap_thermal]
#   type = GapHeatTransfer
#   primary = spacer_facing_upper_plunger
#   secondary = plunger_facing_upper_spacer
#   variable = temperature
#   quadrature = true
#   emissivity_primary = 0.85 #cincotti 2007, table 2
#   emissivity_secondary = 0.85
#   gap_geometry_type = PLATE
#   gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
#   gap_conductivity_function_variable = temperature
#   normal_smoothing_distance = 0.1
# []
# [upper_plunger_diewall_gap_thermal]
#   type = GapHeatTransfer
#   primary = die_wall_facing_upper_plunger
#   secondary = inner_die_wall
#   variable = temperature
#   quadrature = true
#   emissivity_primary = 0.85 #cincotti 2007, table 2
#   emissivity_secondary = 0.85
#   # gap_geometry_type = PLATE # Not for vertical surfaces
#   gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
#   gap_conductivity_function_variable = temperature
#   normal_smoothing_distance = 0.1
# []
# [upper_diewall_spacer_gap_thermal]
#   type = GapHeatTransfer
#   primary = top_die_wall
#   secondary = bottom_upper_spacer
#   variable = temperature
#   quadrature = true
#   emissivity_primary = 0.85 #cincotti 2007, table 2
#   emissivity_secondary = 0.85
#   gap_geometry_type = PLATE
#   gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
#   gap_conductivity_function_variable = temperature
#   normal_smoothing_distance = 0.1
# []
# []

[Materials]
  ## graphite blocks
  [graphite_thermal]
    type = ADGraphiteThermal
    temperature = temperature
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [graphite_density]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = 1.750e3 #in kg/m^3 from Cincotti et al 2007, Table 2, doi:10.1002/aic
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [graphite_electrical_conductivity]
    type = ADGraphiteElectricalConductivity
    temperature = temperature
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []

  #stainless_steel
  [stainless_steel_thermal]
    type = ADStainlessSteelThermal
    temperature = temperature_stainless_steel
    block = stainless_steel
  []
  [stainless_steel_density]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = 8000
    block = stainless_steel
  []
  [stainless_steel_electrical_conductivity]
    type = ADStainlessSteelElectricalConductivity
    temperature = temperature_stainless_steel
    block = stainless_steel
  []

  # harmonic mean of graphite and stainless steel hardness
  [mean_hardness_graphite_stainless_steel]
    type = GraphiteStainlessMeanHardness
    block = 'graphite_spacers upper_plunger lower_plunger die_wall stainless_steel'
  []

  ## yttria powder compact
  [yttria_thermal_conductivity]
    type = ADParsedMaterial
    coupled_variables = 'temperature'
    expression = '3214.46 / (temperature - 147.73)' #in W/(m-K) #Given from Larry's curve fitting, data from Klein and Croft, JAP, v. 38, p. 1603 and UC report "For Computer Heat Conduction Calculations - A compilation of thermal properties data" by A.L. Edwards, UCRL-50589 (1969)
    # coupled_variables = 'thermal_conductivity_aeh'
    # expression = 'thermal_conductivity_aeh' #in W/(m-K) directly, for now
    property_name = 'thermal_conductivity'
    output_properties = thermal_conductivity
    outputs = 'exodus'
    block = powder_compact
  []
  [yttria_specific_heat_capacity]
    type = ADParsedMaterial
    property_name = heat_capacity
    coupled_variables = 'yttria_heat_capacity_volume_avg'
    expression = 'yttria_heat_capacity_volume_avg' #in J/(K-kg)
    # output_properties = yttria_specific_heat_capacity
    # outputs = 'exodus'
    block = powder_compact
  []
  [yttria_density]
    type = ADParsedMaterial
    property_name = 'density'
    coupled_variables = 'yttria_density_volume_avg'
    expression = 'yttria_density_volume_avg'
    # output_properties = yttria_density
    # outputs = 'exodus'
    block = powder_compact
  []
  [electrical_conductivity]
    type = ADParsedMaterial
    #   coupled_variables = 'yttria_sigma_aeh'
    #   expression = 'yttria_sigma_aeh*1.602e8' #converts to units of J/(V^2-m-s)
    property_name = 'electrical_conductivity'
    output_properties = electrical_conductivity
    outputs = 'exodus'
    block = powder_compact
    # type = ADDerivativeParsedMaterial
    # property_name = electrical_conductivity
    coupled_variables = 'temperature'
    constant_names =       'Q_elec  kB            prefactor_solid  initial_porosity'
    constant_expressions = '1.61    8.617343e-5        1.25e-4           0.38'
    expression = '(1-initial_porosity) * prefactor_solid * exp(-Q_elec/kB/temperature) * 1.602e8' # in eV/(nV^2 s nm) per chat with Larry, last term converts to units of J/(V^2-m-s)
  []

  # Material property converter for DiffusionFluxAux object
  [converter]
    type = MaterialADConverter
    ad_props_in = thermal_conductivity
    reg_props_out = nonad_thermal_conductivity
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  automatic_scaling = true
  line_search = 'none'
  # compute_scaling_once = false

  # force running options
  # petsc_options_iname = '-pc_type -snes_linesearch_type -pc_factor_shift_type -pc_factor_shift_amount'
  # petsc_options_value = 'lu       basic                 NONZERO               1e-15'

  # #mechanical contact options
  # petsc_options = '-snes_ksp_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = ' lu       superlu_dist'
  petsc_options = '-snes_converged_reason -ksp_converged_reason'

  nl_forced_its = 1
  nl_rel_tol = 2e-5 #1e-6 #2e-5 for with mechanics #was 1e-10, for temperature only
  nl_abs_tol = 2e-12 #was 1e-12
  nl_max_its = 20
  l_max_its = 50
  dtmin = 1.0e-4

  # num_steps = 1
  end_time = 50 #600 #900 #15 minutes, rule of thumb from Dennis is 10 minutes
  [Quadrature]
    order = FIFTH #required for thermal and mechanical node-face contact
    side_order = SEVENTH
  []
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 0.05
    optimal_iterations = 8
    iteration_window = 2
  []
[]

[Postprocessors]
  [temperature_pp]
    type = AverageNodalVariableValue
    variable = temperature
    block = 'graphite_spacers upper_plunger lower_plunger die_wall powder_compact'
  []
  [temperature_stainless_steel_pp]
    type = AverageNodalVariableValue
    variable = temperature_stainless_steel
    block = stainless_steel
  []
  [graphite_thermal_conductivity]
    type = ElementAverageValue
    variable = thermal_conductivity
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [graphite_electrical_conductivity]
    type = ElementAverageValue
    variable = electrical_conductivity
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [steel_thermal_conductivity]
    type = ElementAverageValue
    variable = thermal_conductivity
    block = stainless_steel
  []
  [steel_electrical_conductivity]
    type = ElementAverageValue
    variable = electrical_conductivity
    block = stainless_steel
  []
  [yttria_thermal_conductivity]
    type = ElementAverageValue
    variable = thermal_conductivity
    block = powder_compact
  []
  [yttria_electrical_conductivity]
    type = ElementAverageValue
    variable = electrical_conductivity
    block = powder_compact
  []
[]

[Outputs]
  csv = true
  exodus = true
  perf_graph = true
  # [ckpt]
  #   type =Checkpoint
  #   time_step_interval = 1
  #   num_files = 2
  # []
[]

(test/tests/interfacekernels/thermal_conductance/thermal_interface_analytic_solution_two_block.i)

[Mesh]
  [line]
    type = GeneratedMeshGenerator
    xmax = 2
    dim = 1
    nx = 400
  []
  [split]
    type = SubdomainBoundingBoxGenerator
    bottom_left = '0 0 0'
    top_right = '1 0 0'
    block_id = 1
    block_name = 'stainless_steel'
    input = line
  []
  [rename_right]
    type = RenameBlockGenerator
    old_block = 0
    new_block = 'graphite'
    input = split
  []
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    input = rename_right
    primary_block = 'stainless_steel'
    paired_block = 'graphite'
    new_boundary = 'ssg_interface'
  []
[]

[Variables]
  [temperature_graphite]
    initial_condition = 300.0
    block = graphite
  []
  [temperature_stainless_steel]
    initial_condition = 300.0
    block = stainless_steel
  []
[]

[AuxVariables]
  [analytic_potential_graphite]
    block = graphite
  []
  [analytic_potential_stainless_steel]
    block = stainless_steel
  []
  [analytic_temperature_graphite]
    block = graphite
  []
  [analytic_temperature_stainless_steel]
    block = stainless_steel
  []
[]

[Kernels]
  [HeatDiff_graphite]
    type = ADHeatConduction
    variable = temperature_graphite
    block = graphite
  []
  [HeatSource_graphite]
    type = ADJouleHeatingSource
    variable = temperature_graphite
    elec = analytic_potential_graphite
    electrical_conductivity = electrical_conductivity
    block = graphite
  []

  [HeatDiff_stainless_steel]
    type = ADHeatConduction
    variable = temperature_stainless_steel
    block = stainless_steel
  []
  [HeatSource_stainless_steel]
    type = ADJouleHeatingSource
    variable = temperature_stainless_steel
    elec = analytic_potential_stainless_steel
    electrical_conductivity = electrical_conductivity
    block = stainless_steel
  []
[]

[AuxKernels]
  [analytic_potential_function_aux_stainless_steel]
    type = FunctionAux
    function = potential_fxn_stainless_steel
    variable = analytic_potential_stainless_steel
    block = stainless_steel
  []
  [analytic_potential_function_aux_graphite]
    type = FunctionAux
    function = potential_fxn_graphite
    variable = analytic_potential_graphite
    block = graphite
  []
  [analytic_temperature_function_aux_stainless_steel]
    type = FunctionAux
    function = temperature_fxn_stainless_steel
    variable = analytic_temperature_stainless_steel
    block = stainless_steel
  []
  [analytic_temperature_function_aux_graphite]
    type = FunctionAux
    function = temperature_fxn_graphite
    variable = analytic_temperature_graphite
    block = graphite
  []
[]

[BCs]
  [thermal_left]
    type = DirichletBC
    variable = temperature_stainless_steel
    boundary = left
    value = 300
  []
  [thermal_right]
    type = DirichletBC
    variable = temperature_graphite
    boundary = right
    value = 300
  []
[]

[InterfaceKernels]
  [thermal_contact_conductance]
    type = ThermalContactCondition
    variable = temperature_stainless_steel
    neighbor_var = temperature_graphite
    primary_potential = analytic_potential_stainless_steel
    secondary_potential = analytic_potential_graphite
    primary_thermal_conductivity = thermal_conductivity
    secondary_thermal_conductivity = thermal_conductivity
    primary_electrical_conductivity = electrical_conductivity
    secondary_electrical_conductivity = electrical_conductivity
    user_electrical_contact_conductance = 75524 # at 300K with 3 kN/m^2 applied pressure
    user_thermal_contact_conductance = 0.242 # also from Cincotti et al
    boundary = ssg_interface
  []
[]

[Materials]
  #graphite
  [heat_conductor_graphite]
    type = ADGenericConstantMaterial
    prop_names = thermal_conductivity
    prop_values = 100
    block = graphite
  []
  [sigma_graphite]
    type = ADGenericConstantMaterial
    prop_names = electrical_conductivity
    prop_values = 73069.2
    block = graphite
  []

  #stainless_steel
  [heat_conductor_stainless_steel]
    type = ADGenericConstantMaterial
    prop_names = thermal_conductivity
    prop_values = 15
    block = stainless_steel
  []
  [sigma_stainless_steel]
    type = ADGenericConstantMaterial
    prop_names = electrical_conductivity
    prop_values = 1.41867e6
    block = stainless_steel
  []
[]

[Functions]
  [potential_fxn_stainless_steel]
    type = ThermalContactPotentialTestFunc
    domain = stainless_steel
  []
  [potential_fxn_graphite]
    type = ThermalContactPotentialTestFunc
    domain = graphite
  []
  [temperature_fxn_stainless_steel]
    type = ThermalContactTemperatureTestFunc
    domain = stainless_steel
  []
  [temperature_fxn_graphite]
    type = ThermalContactTemperatureTestFunc
    domain = graphite
  []
[]

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

[Executioner]
  type = Steady
  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
[]

[Outputs]
  exodus = true
  perf_graph = true
[]

(examples/sps/multiapp/electrothermal_with_phase_field/twoway_initial_prototype/engineering_scale_electrothermal_twoway_prototype.i)

# Approximate Dr. Sinter Geometry
# Configuration: Two stainless steel electrodes (with cooling channels) with two
#                large spacers, and one die-and-plungers assembly.
#                RZ, axisymmetric.
# BCs:
#    Potential:
#       V (top electrode, top surface) --> Neumann condition specifiying potential
#                                          based on applied RMS Current (980 A)
#                                          and cross-sectional electrode area (see
#                                          elec_top BC below). Current turned off
#                                          at 1200 s.
#       V (bottom electrode, bottom surface) = 0 V
#       V (elsewhere) --> natural boundary conditions (no current external to circuit)
#    Temperature:
#       T (top electrode, top surface) = 293 K
#       T (bottom electrode, bottom surface) = 293 K
#       T (vertical right side) --> simple radiative BC into black body at 293 K
#       T (horizonal surfaces) --> GapHeatTransfer
#       T (water channel) --> simple convective BC into "water" with heat transfer coefficient parameter from Cincotti and temperature of 293 K
#       T (elsewhere) --> natural boundary conditions
#    Mechanics (not added yet):
#       Displacement: --> pinned axial centerpoint along interior surface of powder block, outer die wall
#       Pressure: top and bottom surfaces of upper and lower rams, 20.7 MPa
#                 ramped quickly at simulation start, then held constant
# Interface Conditions:
#       V (SS-G upper and G-SS lower interface) --> ElectrostaticContactCondition
#       T (SS-G upper and G-SS lower interface) --> ThermalContactCondition
#       V (graphite-yttria and graphite-graphite) --> GapHeatTransfer
#       T (graphite-yttria and graphite-graphite) --> GapHeatTransfer
#       P (graphite-yttria and graphite-graphite) --> Mechanical Contact
# Initial Conditions:
#       V = default (0 V)
#       T = 873 K (~600C)
# Applied Mechanical Load: needs to be updated to 6.5kN, was 3.5 kN
#
# Reference for graphite, stainless steel: Cincotti et al, DOI 10.1002/aic.11102
# Assorted references for yttria, listed as comments in input file

initial_temperature = 873 #roughly 600C where the pyrometer kicks in
#initial_porosity=0.36 #Maximum random jammed packing, Donev et al (2004) Science Magazine

[Mesh]
  file = ../../../meshes/drsinter_nodalcontact_2d.e
  coord_type = RZ
  construct_side_list_from_node_list = true
  patch_update_strategy = iteration
  patch_size = 20
  # ghosting_patch_size = 5*patch_size
[]

[Problem]
  type = ReferenceResidualProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
[]

[Variables]
  [temperature_stainless_steel]
    initial_condition = 300 #units of K
    block = stainless_steel
  []
  [temperature]
    initial_condition = ${initial_temperature}
    block = 'graphite_spacers lower_plunger upper_plunger die_wall powder_compact'
  []
  [potential_stainless_steel]
    block = stainless_steel
  []
  [electric_potential]
    block = 'graphite_spacers lower_plunger upper_plunger die_wall powder_compact'
  []
[]

[AuxVariables]
  [electric_field_x]
    family = MONOMIAL #prettier pictures with smoother values
    order = FIRST
  []
  [electric_field_y]
    family = MONOMIAL
    order = FIRST
  []
  [current_density]
    family = NEDELEC_ONE
    order = FIRST
  []
  [yttria_sigma_aeh]
    initial_condition = 2.0e-10 #in units eV/((nV)^2-s-nm)
    block = 'powder_compact'
  []
  [microapp_potential]
    #converted to microapp electronVolts units
    block = 'powder_compact'
  []
  [E_x]
    order = FIRST
    family = MONOMIAL
    block = 'powder_compact'
  []
  [E_y]
    order = FIRST
    family = MONOMIAL
    block = 'powder_compact'
  []
  [yttria_current_density_forBC_microapp]
    order = FIRST
    family = MONOMIAL
    block = 'powder_compact'
  []

  [heatflux_graphite_x]
    family = MONOMIAL
    order = CONSTANT
    block = 'graphite_spacers lower_plunger upper_plunger die_wall'
  []
  [heatflux_graphite_y]
    family = MONOMIAL
    order = CONSTANT
    block = 'graphite_spacers lower_plunger upper_plunger die_wall'
  []
  [heatflux_stainless_steel_x]
    family = MONOMIAL
    order = CONSTANT
    block = stainless_steel
  []
  [heatflux_stainless_steel_y]
    family = MONOMIAL
    order = CONSTANT
    block = stainless_steel
  []
  [heatflux_yttria_x]
    family = MONOMIAL
    order = CONSTANT
    block = 'powder_compact'
  []
  [heatflux_yttria_y]
    family = MONOMIAL
    order = CONSTANT
    block = 'powder_compact'
  []
  [yttria_thermal_conductivity_aeh]
    initial_condition = 3.0
  []
  [yttria_heat_capacity_volume_avg]
    initial_condition = 546.914 # at 873K probably? #842.2 at 1600K probably? # at 1500K #568.73 at 1000K #447.281 # at 293K
  []
  [yttria_density_volume_avg]
    initial_condition = 3206.4 ##5010.0*(1-${initial_porosity}) #in kg/m^3
  []
  [heat_transfer_radiation]
  []
  # [thermal_conductivity]
  #   family = MONOMIAL
  #   order = FIRST
  # []
  [Q_from_sub]
    #this will be in eV/m/s, will need unit conversion to J/m^3/s based on phase-field domain size
    order = FIRST
    family = LAGRANGE
  []
[]

[Kernels]
  [HeatDiff_graphite]
    type = ADHeatConduction
    variable = temperature
    thermal_conductivity = thermal_conductivity
    extra_vector_tags = 'ref'
    block = 'graphite_spacers lower_plunger upper_plunger die_wall'
  []
  [HeatTdot_graphite]
    type = ADHeatConductionTimeDerivative
    variable = temperature
    specific_heat = heat_capacity
    density_name = density
    extra_vector_tags = 'ref'
    block = 'graphite_spacers lower_plunger upper_plunger die_wall'
  []
  [JouleHeating_graphite]
    type = ADJouleHeatingSource
    variable = temperature
    elec = electric_potential
    electrical_conductivity = electrical_conductivity
    use_displaced_mesh = true
    extra_vector_tags = 'ref'
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [electric_graphite]
    type = ADMatDiffusion
    variable = electric_potential
    diffusivity = electrical_conductivity
    extra_vector_tags = 'ref'
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []

  [HeatDiff_stainless_steel]
    type = ADHeatConduction
    variable = temperature_stainless_steel
    thermal_conductivity = thermal_conductivity
    extra_vector_tags = 'ref'
    block = stainless_steel
  []
  [HeatTdot_stainless_steel]
    type = ADHeatConductionTimeDerivative
    variable = temperature_stainless_steel
    specific_heat = heat_capacity
    density_name = density
    extra_vector_tags = 'ref'
    block = stainless_steel
  []
  [HeatSource_stainless_steel]
    type = ADJouleHeatingSource
    variable = temperature_stainless_steel
    elec = potential_stainless_steel
    electrical_conductivity = electrical_conductivity
    use_displaced_mesh = true
    extra_vector_tags = 'ref'
    block = stainless_steel
  []
  [electric_stainless_steel]
    type = ADMatDiffusion
    variable = potential_stainless_steel
    diffusivity = electrical_conductivity
    extra_vector_tags = 'ref'
    block = stainless_steel
  []

  [HeatDiff_yttria]
    type = ADHeatConduction
    variable = temperature
    thermal_conductivity = thermal_conductivity #use parsed material property, hope it works
    extra_vector_tags = 'ref'
    block = powder_compact
  []
  [HeatTdot_yttria]
    type = ADHeatConductionTimeDerivative
    variable = temperature
    specific_heat = heat_capacity #use parsed material property
    density_name = density
    extra_vector_tags = 'ref'
    block = powder_compact
  []
  [heat_source]
    type = MaskedBodyForce
    mask = Q_SI
    variable = temperature
  []
  # [JouleHeating_yttria]
  #   type = ADJouleHeatingSource
  #   variable = temperature
  #   elec = electric_potential
  #   electrical_conductivity = electrical_conductivity
  #   use_displaced_mesh = true
  #   extra_vector_tags = 'ref'
  #   block = powder_compact
  # []
  [electric_yttria]
    type = ADMatDiffusion
    variable = electric_potential
    diffusivity = electrical_conductivity
    extra_vector_tags = 'ref'
    block = powder_compact
  []
[]

[AuxKernels]
  [electrostatic_calculation_x_graphite_yttria]
    type = PotentialToFieldAux
    gradient_variable = electric_potential
    variable = electric_field_x
    sign = negative
    component = x
    block = 'powder_compact graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [electrostatic_calculation_y_graphite_yttria]
    type = PotentialToFieldAux
    gradient_variable = electric_potential
    variable = electric_field_y
    sign = negative
    component = y
    block = 'powder_compact graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [electrostatic_calculation_x_stainless_steel]
    type = PotentialToFieldAux
    gradient_variable = potential_stainless_steel
    variable = electric_field_x
    sign = negative
    component = x
    block = stainless_steel
  []
  [electrostatic_calculation_y_stainless_steel]
    type = PotentialToFieldAux
    gradient_variable = potential_stainless_steel
    variable = electric_field_y
    sign = negative
    component = y
    block = stainless_steel
  []
  [current_density_graphite_yttria]
    type = ADCurrentDensity
    variable = current_density
    potential = electric_potential
    block = 'powder_compact graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [current_density_stainless_steel]
    type = ADCurrentDensity
    variable = current_density
    potential = potential_stainless_steel
    block = stainless_steel
  []
  [microapp_potential]
    type = ParsedAux
    variable = microapp_potential
    coupled_variables = electric_potential
    expression = 'electric_potential*1e9' #convert from V to nV
    block = 'powder_compact'
  []
  [E_x]
    type = VariableGradientComponent
    variable = E_x
    gradient_variable = electric_potential
    component = x
    block = 'powder_compact'
  []
  [E_y]
    type = VariableGradientComponent
    variable = E_y
    gradient_variable = electric_potential
    component = y
    block = 'powder_compact'
  []
  [yttria_current_density_forBC_microapp]
    type = ParsedAux
    variable = yttria_current_density_forBC_microapp
    coupled_variables = 'electrical_conductivity E_y'
    expression = '-1.0*electrical_conductivity*E_y'
    block = 'powder_compact'
  []

  [heat_flux_graphite_x]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature
    variable = heatflux_graphite_x
    component = x
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [heat_flux_graphite_y]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature
    variable = heatflux_graphite_y
    component = y
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [heat_flux_stainless_x]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature_stainless_steel
    variable = heatflux_stainless_steel_x
    block = stainless_steel
    component = x
  []
  [heat_flux_stainless_y]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature_stainless_steel
    variable = heatflux_stainless_steel_y
    block = stainless_steel
    component = y
  []
  [heat_flux_yttria_x]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature
    variable = heatflux_yttria_x
    component = x
    block = 'powder_compact'
  []
  [heat_flux_yttria_y]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature
    variable = heatflux_yttria_y
    component = y
    block = 'powder_compact'
  []

  [heat_transfer_radiation_graphite]
    type = ParsedAux
    variable = heat_transfer_radiation
    boundary = 'outer_radiative_spacers outer_die_wall radiative_upper_plunger radiative_lower_plunger'
    coupled_variables = 'temperature'
    constant_names = 'boltzmann epsilon temperature_farfield' #published emissivity for graphite is 0.85
    constant_expressions = '5.67e-8 0.85 293.0' #roughly room temperature, which is probably too cold
    expression = '-boltzmann*epsilon*(temperature^4-temperature_farfield^4)'
  []
  [heat_transfer_radiation_stainless_steel]
    type = ParsedAux
    variable = heat_transfer_radiation
    boundary = 'outer_radiative_stainless_steel'
    coupled_variables = 'temperature_stainless_steel'
    constant_names = 'boltzmann epsilon temperature_farfield' #published emissivity for graphite is 0.85
    constant_expressions = '5.67e-8 0.4 293.0' #roughly room temperature, which is probably too cold
    expression = '-boltzmann*epsilon*(temperature_stainless_steel^4-temperature_farfield^4)'
  []
  [thermal_conductivity_graphite]
    type = ADMaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [thermal_conductivity_steel]
    type = ADMaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = 'stainless_steel'
  []
  [thermal_conductivity_yttria]
    type = ADMaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = 'powder_compact'
  []
  [electrical_conductivity_graphite]
    type = ADMaterialRealAux
    variable = electrical_conductivity
    property = electrical_conductivity
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [electrical_conductivity_steel]
    type = ADMaterialRealAux
    variable = electrical_conductivity
    property = electrical_conductivity
    block = 'stainless_steel'
  []
  [electrical_conductivity_yttria]
    type = ADMaterialRealAux
    variable = electrical_conductivity
    property = electrical_conductivity
    block = 'powder_compact'
  []
[]

[BCs]
  [external_surface_stainless]
    type = CoupledVarNeumannBC
    boundary = 'outer_radiative_stainless_steel'
    variable = temperature_stainless_steel
    v = heat_transfer_radiation
  []
  [external_surface_graphite]
    type = CoupledVarNeumannBC
    boundary = 'outer_radiative_spacers outer_die_wall radiative_upper_plunger radiative_lower_plunger'
    variable = temperature
    v = heat_transfer_radiation
  []
  [water_channel]
    type = CoupledConvectiveHeatFluxBC
    boundary = water_channel
    variable = temperature_stainless_steel
    T_infinity = 300 #units of K
    htc = 4725
  []
  [temperature_ram_extremes]
    type = ADDirichletBC
    variable = temperature_stainless_steel
    boundary = 'top_upper_steel_ram bottom_lower_steel_ram'
    value = 293
  []
  [electric_potential_top]
    type = ADFunctionNeumannBC
    variable = potential_stainless_steel
    boundary = top_upper_steel_ram
    ## This will need to be updated to match the Dr. Sinter geometry
    function = 'if(t < 31, (980 / (pi * 0.00155))*((0.141301/4.3625)*t), if(t > 600, 0, 980 / (pi * 0.00155)))' # RMS Current / Cross-sectional Area. Ramping for t < 31s approximately reflects Cincotti et al (DOI: 10.1002/aic.11102) Figure 21(b)
  []
  [electric_potential_bottom]
    type = ADDirichletBC
    variable = potential_stainless_steel
    boundary = bottom_lower_steel_ram
    value = 0
  []
[]

[Functions]
  [graphite_thermal_conductivity_fcn]
    type = ParsedFunction
    expression = '-4.418e-12*t^4+2.904e-8*t^3-4.688e-5*t^2-0.0316*t+119.659'
  []
  # [yttria_thermal_conductivity_fcn] #from the multiapp
  #   type = ParsedFunction
  #   expression = '3214.46/(t-147.73)'
  # []
  [harmonic_mean_thermal_conductivity]
    type = ParsedFunction
    expression = '2*(-4.418e-12*t^4+2.904e-8*t^3-4.688e-5*t^2-0.0316*t+119.659)*(3214.46/(t-147.73))/((-4.418e-12*t^4+2.904e-8*t^3-4.688e-5*t^2-0.0316*t+119.659)+(3214.46/(t-147.73)))'
    # symbol_names = 'k_graphite k_yttria'
    # symbol_values = 'graphite_thermal_conductivity_fcn yttria_thermal_conductivity_fcn'
  []

  [graphite_electrical_conductivity_fcn]
    type = ParsedFunction
    expression = '1.0/(-2.705e-15*t^3+1.263e-11*t^2-1.836e-8*t+1.813e-5)'
  []
  # [electrical_conductivity_fcn]
  #   type = ParsedFunction
  #   # symbol_names = porosity
  #   # symbol_values = initial_porosity
  #   expression = '(1-0.62)*2.0025e4*exp(-1.61/8.617343e-5/t)'
  # []
  [harmonic_mean_electrical_conductivity]
    type = ParsedFunction
    expression = '2*(1.0/(-2.705e-15*t^3+1.263e-11*t^2-1.836e-8*t+1.813e-5))*((1)*2.0025e4*exp(-1.61/8.617343e-5/t))/((1.0/(-2.705e-15*t^3+1.263e-11*t^2-1.836e-8*t+1.813e-5))+((1)*2.0025e4*exp(-1.61/8.617343e-5/t)))'
    # symbol_names = 'k_graphite k_yttria'
    # symbol_values = 'graphite_thermal_conductivity_fcn yttria_thermal_conductivity_fcn'
  []
  [mechanical_pressure_func]
    type = ParsedFunction
    symbol_names = 'radius coolant_radius force'
    symbol_values = '0.04 7.071e-3 6.5' # 'm kN'
    expression = 'force * 1e3 / (pi * (radius^2 - coolant_radius^2))' # (N / m^2)
  []
[]

[InterfaceKernels]
  [electric_contact_conductance_ssg]
    type = ElectrostaticContactCondition
    variable = potential_stainless_steel
    neighbor_var = electric_potential
    boundary = 'upper_ram_spacer_interface lower_ram_spacer_interface'
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
  []
  [thermal_contact_conductance_calculated_ssg]
    type = ThermalContactCondition
    variable = temperature_stainless_steel
    neighbor_var = temperature
    primary_potential = potential_stainless_steel
    secondary_potential = electric_potential
    primary_thermal_conductivity = thermal_conductivity
    secondary_thermal_conductivity = thermal_conductivity
    primary_electrical_conductivity = electrical_conductivity
    secondary_electrical_conductivity = electrical_conductivity
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
    boundary = 'upper_ram_spacer_interface lower_ram_spacer_interface'
  []
  [electric_contact_conductance_gss]
    type = ElectrostaticContactCondition
    variable = electric_potential
    neighbor_var = potential_stainless_steel
    boundary = 'upper_ram_spacer_interface lower_ram_spacer_interface'
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
  []
  [thermal_contact_conductance_calculated_gss]
    type = ThermalContactCondition
    variable = temperature
    neighbor_var = temperature_stainless_steel
    primary_potential = electric_potential
    secondary_potential = potential_stainless_steel
    primary_thermal_conductivity = thermal_conductivity
    secondary_thermal_conductivity = thermal_conductivity
    primary_electrical_conductivity = electrical_conductivity
    secondary_electrical_conductivity = electrical_conductivity
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
    boundary = 'upper_ram_spacer_interface lower_ram_spacer_interface'
  []
[]

[ThermalContact]
  [upper_plunger_powder_electric]
    type = GapHeatTransfer
    primary = bottom_upper_plunger
    secondary = top_powder_compact
    variable = electric_potential
    quadrature = true
    emissivity_primary = 0.0 #not applicable for electric potential
    emissivity_secondary = 0.0 #not applicable for electric potential
    gap_geometry_type = PLATE
    # gap_conductivity = 1e-23 #rough harmonic mean
    gap_conductivity_function = 'harmonic_mean_electrical_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [powder_bottom_plunger_electric]
    type = GapHeatTransfer
    primary = top_lower_plunger
    secondary = bottom_powder_compact
    variable = electric_potential
    quadrature = true
    emissivity_primary = 0.0 #not applicable for electric potential
    emissivity_secondary = 0.0 #not applicable for electric potential
    gap_geometry_type = PLATE
    # gap_conductivity = 1e-23 #rough harmonic mean
    gap_conductivity_function = 'harmonic_mean_electrical_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [powder_die_electric]
    type = GapHeatTransfer
    primary = inner_die_wall
    secondary = outer_powder_compact
    variable = electric_potential
    quadrature = true
    emissivity_primary = 0.0 #0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.0 #0.85
    # gap_geometry_type = PLATE # Not for vertical surfaces
    # gap_conductivity = 1e-23 #rough harmonic mean
    gap_conductivity_function = 'harmonic_mean_electrical_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [combined_plungers_die_electric]
    type = GapHeatTransfer
    primary = inner_die_wall
    secondary = combined_outer_plungers
    variable = electric_potential
    quadrature = true
    emissivity_primary = 0.0 #not applicable for electric potential
    emissivity_secondary = 0.0 #not applicable for electric potential
    # gap_geometry_type = PLATE # Not for vertical surfaces
    # gap_conductivity = 8.5e4 #graphite, at 500K
    gap_conductivity_function = 'graphite_electrical_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []

##Thermal Contact between gapped graphite die components
  [upper_plunger_spacer_gap_thermal]
    type = GapHeatTransfer
    primary = spacer_facing_upper_plunger
    secondary = plunger_facing_upper_spacer
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    gap_geometry_type = PLATE
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [upper_diewall_spacer_gap_thermal]
    type = GapHeatTransfer
    primary = top_die_wall
    secondary = bottom_upper_spacer
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    gap_geometry_type = PLATE
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [lower_plunger_spacer_gap_thermal]
    type = GapHeatTransfer
    primary = spacer_facing_lower_plunger
    secondary = plunger_facing_lower_spacer
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    gap_geometry_type = PLATE
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [lower_diewall_spacer_gap_thermal]
    type = GapHeatTransfer
    primary = bottom_die_wall
    secondary = top_lower_spacer
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    gap_geometry_type = PLATE
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [upper_plunger_diewall_gap_thermal]
    type = GapHeatTransfer
    primary = inner_die_wall ### paired temperature doesn't show on inner die wall, but temperature profile looks reasonable
    secondary = die_wall_facing_upper_plunger
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    # gap_geometry_type = PLATE # Not for vertical surfaces
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [lower_plunger_diewall_gap_thermal]
    type = GapHeatTransfer
    primary = inner_die_wall ### paired temperature doesn't show on inner die wall, but temperature profile looks reasonable
    secondary = die_wall_facing_lower_plunger
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    # gap_geometry_type = PLATE # Not for vertical surfaces
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []

## Thermal Contact between touching components of powder and die
  [upper_plunger_powder_thermal]
    type = GapHeatTransfer
    primary = bottom_upper_plunger
    secondary = top_powder_compact
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.3 #estimated from McMahon and Wilder, High Temperature Spectral Emissivity of Yttrium, Samarium, Gadolinium, Ebrium and Lutetium Oxides (1963) Atomic Energy Commission, IS-578, Figure 12
    gap_geometry_type = PLATE
    gap_conductivity_function = 'harmonic_mean_thermal_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [powder_bottom_plunger_thermal]
    type = GapHeatTransfer
    primary = top_lower_plunger
    secondary = bottom_powder_compact #expect more heat transfer from the die to the powder
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.3 #estimated from McMahon and Wilder, High Temperature Spectral Emissivity of Yttrium, Samarium, Gadolinium, Ebrium and Lutetium Oxides (1963) Atomic Energy Commission, IS-578, Figure 12
    gap_geometry_type = PLATE
    gap_conductivity_function = 'harmonic_mean_thermal_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [powder_die_thermal]
    type = GapHeatTransfer
    primary = inner_die_wall
    secondary = outer_powder_compact
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.3 #estimated from McMahon and Wilder, High Temperature Spectral Emissivity of Yttrium, Samarium, Gadolinium, Ebrium and Lutetium Oxides (1963) Atomic Energy Commission, IS-578, Figure 12
    # gap_geometry_type = PLATE # Not for vertical surfaces
    gap_conductivity_function = 'harmonic_mean_thermal_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [combined_plungers_die_thermal]
    type = GapHeatTransfer
    primary = inner_die_wall
    secondary = combined_outer_plungers
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    # gap_geometry_type = PLATE # Not for vertical surfaces
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
[]

[Materials]
  ## graphite blocks
  [graphite_thermal]
    type = ADGraphiteThermal
    temperature = temperature
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [graphite_density]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = 1.750e3 #in kg/m^3 from Cincotti et al 2007, Table 2, doi:10.1002/aic
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [graphite_electrical_conductivity]
    type = ADGraphiteElectricalConductivity
    temperature = temperature
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []

  #stainless_steel
  [stainless_steel_thermal]
    type = ADStainlessSteelThermal
    temperature = temperature_stainless_steel
    block = stainless_steel
  []
  [stainless_steel_density]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = 8000
    block = stainless_steel
  []
  [stainless_steel_electrical_conductivity]
    type = ADStainlessSteelElectricalConductivity
    temperature = temperature_stainless_steel
    block = stainless_steel
  []

  # harmonic mean of graphite and stainless steel hardness
  [mean_hardness_graphite_stainless_steel]
    type = GraphiteStainlessMeanHardness
    block = 'graphite_spacers upper_plunger lower_plunger die_wall stainless_steel'
  []

  ## yttria powder compact
  [yttria_thermal_conductivity]
    type = ADParsedMaterial
    # coupled_variables = 'temperature'
    # expression = '3214.46 / (temperature - 147.73)' #in W/(m-K) #Given from Larry's curve fitting, data from Klein and Croft, JAP, v. 38, p. 1603 and UC report "For Computer Heat Conduction Calculations - A compilation of thermal properties data" by A.L. Edwards, UCRL-50589 (1969)
    coupled_variables = 'yttria_thermal_conductivity_aeh'
    expression = 'yttria_thermal_conductivity_aeh' #in W/(m-K) directly, for now
    property_name = 'thermal_conductivity'
    output_properties = thermal_conductivity
    outputs = 'exodus'
    block = powder_compact
  []
  [yttria_specific_heat_capacity]
    type = ADParsedMaterial
    property_name = heat_capacity
    coupled_variables = 'yttria_heat_capacity_volume_avg'
    expression = 'yttria_heat_capacity_volume_avg' #in J/(K-kg)
    # output_properties = yttria_specific_heat_capacity
    # outputs = 'exodus'
    block = powder_compact
  []
  [yttria_density]
    type = ADParsedMaterial
    property_name = 'density'
    coupled_variables = 'yttria_density_volume_avg'
    expression = 'yttria_density_volume_avg'
    # output_properties = yttria_density
    # outputs = 'exodus'
    block = powder_compact
  []
  [electrical_conductivity]
    type = ADParsedMaterial
    coupled_variables = 'yttria_sigma_aeh'
    expression = 'yttria_sigma_aeh*1.602e-10' #converts from eV/(V^2 s nm) to units of J/(V^2-m-s)
    property_name = 'electrical_conductivity'
    output_properties = electrical_conductivity
    outputs = 'exodus'
    block = powder_compact
    # type = ADDerivativeParsedMaterial
    # property_name = electrical_conductivity
    # coupled_variables = 'temperature'
    # constant_names =       'Q_elec  kB            prefactor_solid  initial_porosity'
    # constant_expressions = '1.61    8.617343e-5        1.25e-4           0.38'
    # expression = '(1-initial_porosity) * prefactor_solid * exp(-Q_elec/kB/temperature) * 1.602e-10' # in eV/(V^2 s nm) per chat with Larry, last term converts to units of J/(V^2-m-s)
  []

  [Q_SI]
    type = ParsedMaterial
    property_name = Q_SI
    coupled_variables = 'Q_from_sub'
    expression = 'Q_from_sub / 80 / 40 * 0.1602' #divide by domain size and unit conversion to go from eV/S to J/m^3/s
    outputs = 'exodus'
  []

  # Material property converter for DiffusionFluxAux object
  [converter]
    type = MaterialADConverter
    ad_props_in = thermal_conductivity
    reg_props_out = nonad_thermal_conductivity
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  automatic_scaling = true
  line_search = 'none'
  # compute_scaling_once = false

  # force running options
  # petsc_options_iname = '-pc_type -snes_linesearch_type -pc_factor_shift_type -pc_factor_shift_amount'
  # petsc_options_value = 'lu       basic                 NONZERO               1e-15'

  # #mechanical contact options
  # petsc_options = '-snes_ksp_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = ' lu       superlu_dist'
  petsc_options = '-snes_converged_reason -ksp_converged_reason'

  nl_forced_its = 1
  nl_rel_tol = 2e-5 #1e-6 #2e-5 for with mechanics #was 1e-10, for temperature only
  nl_abs_tol = 2e-12 #was 1e-12
  nl_max_its = 20
  l_max_its = 50
  dtmin = 1.0e-4

  end_time = 900 #600 #900 #15 minutes, rule of thumb from Dennis is 10 minutes
  [Quadrature]
    order = FIFTH #required for thermal and mechanical node-face contact
    side_order = SEVENTH
  []
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 0.05
    optimal_iterations = 8
    iteration_window = 2
  []
[]

[Postprocessors]
  [temperature_pp]
    type = AverageNodalVariableValue
    variable = temperature
    block = 'graphite_spacers upper_plunger lower_plunger die_wall powder_compact'
  []
  [temperature_stainless_steel_pp]
    type = AverageNodalVariableValue
    variable = temperature_stainless_steel
    block = stainless_steel
  []
  [graphite_thermal_conductivity]
    type = ElementAverageValue
    variable = thermal_conductivity
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [graphite_electrical_conductivity]
    type = ElementAverageValue
    variable = electrical_conductivity
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [steel_thermal_conductivity]
    type = ElementAverageValue
    variable = thermal_conductivity
    block = stainless_steel
  []
  [steel_electrical_conductivity]
    type = ElementAverageValue
    variable = electrical_conductivity
    block = stainless_steel
  []
  [yttria_thermal_conductivity]
    type = ElementAverageValue
    variable = thermal_conductivity
    block = powder_compact
  []
  [yttria_electrical_conductivity]
    type = ElementAverageValue
    variable = electrical_conductivity
    block = powder_compact
  []
[]

[MultiApps]
  [micro]
    type = TransientMultiApp
    # type = CentroidMultiApp # lauches one in the middle of each element so don't need to give positions
    #can specify the number of procs
    max_procs_per_app = 1 #paolo recommends starting here
    app_type = MalamuteApp
    positions = '0.00125 0.034 0' #roughly the center of element 117 in this mesh
    # positions = '0.0074 0.0058 0' #roughly the center of element 368 in this mesh
    input_files = micro_yttria_thermoelectric_twoway.i
    sub_cycling = true
    execute_on = TIMESTEP_BEGIN #the default
  []
[]

[Transfers]
  [keff_from_sub]
    type = MultiAppPostprocessorInterpolationTransfer
    from_multi_app = micro
    variable = yttria_thermal_conductivity_aeh
    power = 2 #2 is the default value, tutorial uses 1
    postprocessor = k_AEH_average
  []
  [sigma_aeh_eff_from_sub]
    type = MultiAppPostprocessorInterpolationTransfer
    from_multi_app = micro
    variable = yttria_sigma_aeh
    power = 2 #2 is the default value, tutorial uses 1
    postprocessor = sigma_y_AEH
  []
  [Q_from_sub]
    type = MultiAppPostprocessorInterpolationTransfer
    from_multi_app = micro
    variable = Q_from_sub #This is the integrated heat produced in the phase-field simulation in eV/m/s
    power = 2 #2 is the default value, tutorial uses 1
    postprocessor = Q_joule_total
  []
  [temperature_to_sub]
    type = MultiAppVariableValueSampleTransfer
    to_multi_app = micro
    source_variable = temperature
    variable = T
  []
  [temperature_to_sub_postproc]
    type = MultiAppVariableValueSamplePostprocessorTransfer
    to_multi_app = micro
    source_variable = temperature
    postprocessor = T_postproc
  []
  [potential_to_sub_postproc]
    type = MultiAppVariableValueSamplePostprocessorTransfer
    to_multi_app = micro
    source_variable = electric_potential
    postprocessor = V_postproc
  []
  [micro_field_pp_to_sub]
    type = MultiAppVariableValueSamplePostprocessorTransfer
    to_multi_app = micro
    source_variable = E_y
    postprocessor = Ey_in
  []
[]

[Outputs]
  csv = true
  exodus = true
  perf_graph = true
  # [ckpt]
  #   type =Checkpoint
  #   time_step_interval = 1
  #   num_files = 2
  # []
[]

(examples/sps/multiapp/electrothermal_with_phase_field/oneway_controls/engineering_scale_electrothermal_oneway_controls.i)

# Approximate Dr. Sinter Geometry
# Configuration: Two stainless steel electrodes (with cooling channels) with two
#                large spacers, and one die-and-plungers assembly.
#                RZ, axisymmetric.
# BCs:
#    Potential:
#       V (top electrode, top surface) --> Neumann condition specifiying potential
#                                          based on applied RMS Current (980 A)
#                                          and cross-sectional electrode area (see
#                                          elec_top BC below). Current turned off
#                                          at 1200 s.
#       V (bottom electrode, bottom surface) = 0 V
#       V (elsewhere) --> natural boundary conditions (no current external to circuit)
#    Temperature:
#       T (top electrode, top surface) = 293 K
#       T (bottom electrode, bottom surface) = 293 K
#       T (vertical right side) --> simple radiative BC into black body at 293 K
#       T (horizonal surfaces) --> GapHeatTransfer
#       T (water channel) --> simple convective BC into "water" with heat transfer coefficient parameter from Cincotti and temperature of 293 K
#       T (elsewhere) --> natural boundary conditions
#    Mechanics (not added yet):
#       Displacement: --> pinned axial centerpoint along interior surface of powder block, outer die wall
#       Pressure: top and bottom surfaces of upper and lower rams, 20.7 MPa
#                 ramped quickly at simulation start, then held constant
# Interface Conditions:
#       V (SS-G upper and G-SS lower interface) --> ElectrostaticContactCondition
#       T (SS-G upper and G-SS lower interface) --> ThermalContactCondition
#       V (graphite-yttria and graphite-graphite) --> GapHeatTransfer
#       T (graphite-yttria and graphite-graphite) --> GapHeatTransfer
#       P (graphite-yttria and graphite-graphite) --> Mechanical Contact
# Initial Conditions:
#       V = default (0 V)
#       T = 873 K (~600C)
# Applied Mechanical Load: needs to be updated to 6.5kN, was 3.5 kN
#
# Reference for graphite, stainless steel: Cincotti et al, DOI 10.1002/aic.11102
# Assorted references for yttria, listed as comments in input file

initial_temperature = 873 #roughly 600C where the pyrometer kicks in
#initial_porosity=0.36 #Maximum random jammed packing, Donev et al (2004) Science Magazine

[Mesh]
  file = ../../../meshes/drsinter_nodalcontact_2d.e
  coord_type = RZ
  construct_side_list_from_node_list = true
  patch_update_strategy = iteration
  patch_size = 20
  # ghosting_patch_size = 5*patch_size
[]

[Problem]
  type = ReferenceResidualProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
[]

[Variables]
  [temperature_stainless_steel]
    initial_condition = 300 #units of K
    block = stainless_steel
  []
  [temperature]
    initial_condition = ${initial_temperature}
    block = 'graphite_spacers lower_plunger upper_plunger die_wall powder_compact'
  []
  [potential_stainless_steel]
    block = stainless_steel
  []
  [electric_potential]
    block = 'graphite_spacers lower_plunger upper_plunger die_wall powder_compact'
  []
[]

[AuxVariables]
  [electric_field_x]
    family = MONOMIAL #prettier pictures with smoother values
    order = FIRST
  []
  [electric_field_y]
    family = MONOMIAL
    order = FIRST
  []
  [current_density]
    family = NEDELEC_ONE
    order = FIRST
  []
  [yttria_sigma_aeh]
    initial_condition = 2.0e-10 #in units eV/((nV)^2-s-nm)
    block = 'powder_compact'
  []
  [microapp_potential]
    #converted to microapp electronVolts units
    block = 'powder_compact'
  []
  [E_x]
    order = FIRST
    family = MONOMIAL
    block = 'powder_compact'
  []
  [E_y]
    order = FIRST
    family = MONOMIAL
    block = 'powder_compact'
  []
  [yttria_current_density_forBC_microapp]
    order = FIRST
    family = MONOMIAL
    block = 'powder_compact'
  []

  [heatflux_graphite_x]
    family = MONOMIAL
    order = CONSTANT
    block = 'graphite_spacers lower_plunger upper_plunger die_wall'
  []
  [heatflux_graphite_y]
    family = MONOMIAL
    order = CONSTANT
    block = 'graphite_spacers lower_plunger upper_plunger die_wall'
  []
  [heatflux_stainless_steel_x]
    family = MONOMIAL
    order = CONSTANT
    block = stainless_steel
  []
  [heatflux_stainless_steel_y]
    family = MONOMIAL
    order = CONSTANT
    block = stainless_steel
  []
  [heatflux_yttria_x]
    family = MONOMIAL
    order = CONSTANT
    block = 'powder_compact'
  []
  [heatflux_yttria_y]
    family = MONOMIAL
    order = CONSTANT
    block = 'powder_compact'
  []
  [yttria_thermal_conductivity_aeh]
    initial_condition = 0.4
  []
  [yttria_heat_capacity_volume_avg]
    initial_condition = 546.914 # at 873K probably? #842.2 at 1600K probably? # at 1500K #568.73 at 1000K #447.281 # at 293K
  []
  [yttria_density_volume_avg]
    initial_condition = 3106.2 ##5010.0*(1-${initial_porosity}) #in kg/m^3
  []
  [heat_transfer_radiation]
  []
  # [thermal_conductivity]
  #   family = MONOMIAL
  #   order = FIRST
  # []
[]

[Kernels]
  [HeatDiff_graphite]
    type = ADHeatConduction
    variable = temperature
    thermal_conductivity = thermal_conductivity
    extra_vector_tags = 'ref'
    block = 'graphite_spacers lower_plunger upper_plunger die_wall'
  []
  [HeatTdot_graphite]
    type = ADHeatConductionTimeDerivative
    variable = temperature
    specific_heat = heat_capacity
    density_name = density
    extra_vector_tags = 'ref'
    block = 'graphite_spacers lower_plunger upper_plunger die_wall'
  []
  [JouleHeating_graphite]
    type = ADJouleHeatingSource
    variable = temperature
    elec = electric_potential
    electrical_conductivity = electrical_conductivity
    use_displaced_mesh = true
    extra_vector_tags = 'ref'
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [electric_graphite]
    type = ADMatDiffusion
    variable = electric_potential
    diffusivity = electrical_conductivity
    extra_vector_tags = 'ref'
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []

  [HeatDiff_stainless_steel]
    type = ADHeatConduction
    variable = temperature_stainless_steel
    thermal_conductivity = thermal_conductivity
    extra_vector_tags = 'ref'
    block = stainless_steel
  []
  [HeatTdot_stainless_steel]
    type = ADHeatConductionTimeDerivative
    variable = temperature_stainless_steel
    specific_heat = heat_capacity
    density_name = density
    extra_vector_tags = 'ref'
    block = stainless_steel
  []
  [HeatSource_stainless_steel]
    type = ADJouleHeatingSource
    variable = temperature_stainless_steel
    elec = potential_stainless_steel
    electrical_conductivity = electrical_conductivity
    use_displaced_mesh = true
    extra_vector_tags = 'ref'
    block = stainless_steel
  []
  [electric_stainless_steel]
    type = ADMatDiffusion
    variable = potential_stainless_steel
    diffusivity = electrical_conductivity
    extra_vector_tags = 'ref'
    block = stainless_steel
  []

  [HeatDiff_yttria]
    type = ADHeatConduction
    variable = temperature
    thermal_conductivity = thermal_conductivity #use parsed material property, hope it works
    extra_vector_tags = 'ref'
    block = powder_compact
  []
  [HeatTdot_yttria]
    type = ADHeatConductionTimeDerivative
    variable = temperature
    specific_heat = heat_capacity #use parsed material property
    density_name = density
    extra_vector_tags = 'ref'
    block = powder_compact
  []
  [JouleHeating_yttria]
    type = ADJouleHeatingSource
    variable = temperature
    elec = electric_potential
    electrical_conductivity = electrical_conductivity
    use_displaced_mesh = true
    extra_vector_tags = 'ref'
    block = powder_compact
  []
  [electric_yttria]
    type = ADMatDiffusion
    variable = electric_potential
    diffusivity = electrical_conductivity
    extra_vector_tags = 'ref'
    block = powder_compact
  []
[]

[AuxKernels]
  [electrostatic_calculation_x_graphite_yttria]
    type = PotentialToFieldAux
    gradient_variable = electric_potential
    variable = electric_field_x
    sign = negative
    component = x
    block = 'powder_compact graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [electrostatic_calculation_y_graphite_yttria]
    type = PotentialToFieldAux
    gradient_variable = electric_potential
    variable = electric_field_y
    sign = negative
    component = y
    block = 'powder_compact graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [electrostatic_calculation_x_stainless_steel]
    type = PotentialToFieldAux
    gradient_variable = potential_stainless_steel
    variable = electric_field_x
    sign = negative
    component = x
    block = stainless_steel
  []
  [electrostatic_calculation_y_stainless_steel]
    type = PotentialToFieldAux
    gradient_variable = potential_stainless_steel
    variable = electric_field_y
    sign = negative
    component = y
    block = stainless_steel
  []
  [current_density_graphite_yttria]
    type = ADCurrentDensity
    variable = current_density
    potential = electric_potential
    block = 'powder_compact graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [current_density_stainless_steel]
    type = ADCurrentDensity
    variable = current_density
    potential = potential_stainless_steel
    block = stainless_steel
  []
  [microapp_potential]
    type = ParsedAux
    variable = microapp_potential
    coupled_variables = electric_potential
    expression = 'electric_potential*1e9' #convert from V to nV
    block = 'powder_compact'
  []
  [E_x]
    type = VariableGradientComponent
    variable = E_x
    gradient_variable = electric_potential
    component = x
    block = 'powder_compact'
  []
  [E_y]
    type = VariableGradientComponent
    variable = E_y
    gradient_variable = electric_potential
    component = y
    block = 'powder_compact'
  []
  [yttria_current_density_forBC_microapp]
    type = ParsedAux
    variable = yttria_current_density_forBC_microapp
    coupled_variables = 'electrical_conductivity E_y'
    expression = '-1.0*electrical_conductivity*E_y'
    block = 'powder_compact'
  []

  [heat_flux_graphite_x]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature
    variable = heatflux_graphite_x
    component = x
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [heat_flux_graphite_y]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature
    variable = heatflux_graphite_y
    component = y
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [heat_flux_stainless_x]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature_stainless_steel
    variable = heatflux_stainless_steel_x
    block = stainless_steel
    component = x
  []
  [heat_flux_stainless_y]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature_stainless_steel
    variable = heatflux_stainless_steel_y
    block = stainless_steel
    component = y
  []
  [heat_flux_yttria_x]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature
    variable = heatflux_yttria_x
    component = x
    block = 'powder_compact'
  []
  [heat_flux_yttria_y]
    type = DiffusionFluxAux
    diffusivity = nonad_thermal_conductivity
    diffusion_variable = temperature
    variable = heatflux_yttria_y
    component = y
    block = 'powder_compact'
  []

  [heat_transfer_radiation_graphite]
    type = ParsedAux
    variable = heat_transfer_radiation
    boundary = 'outer_radiative_spacers outer_die_wall radiative_upper_plunger radiative_lower_plunger'
    coupled_variables = 'temperature'
    constant_names = 'boltzmann epsilon temperature_farfield' #published emissivity for graphite is 0.85
    constant_expressions = '5.67e-8 0.85 293.0' #roughly room temperature, which is probably too cold
    expression = '-boltzmann*epsilon*(temperature^4-temperature_farfield^4)'
  []
  [heat_transfer_radiation_stainless_steel]
    type = ParsedAux
    variable = heat_transfer_radiation
    boundary = 'outer_radiative_stainless_steel'
    coupled_variables = 'temperature_stainless_steel'
    constant_names = 'boltzmann epsilon temperature_farfield' #published emissivity for graphite is 0.85
    constant_expressions = '5.67e-8 0.4 293.0' #roughly room temperature, which is probably too cold
    expression = '-boltzmann*epsilon*(temperature_stainless_steel^4-temperature_farfield^4)'
  []
  [thermal_conductivity_graphite]
    type = ADMaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [thermal_conductivity_steel]
    type = ADMaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = 'stainless_steel'
  []
  [thermal_conductivity_yttria]
    type = ADMaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = 'powder_compact'
  []
  [electrical_conductivity_graphite]
    type = ADMaterialRealAux
    variable = electrical_conductivity
    property = electrical_conductivity
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [electrical_conductivity_steel]
    type = ADMaterialRealAux
    variable = electrical_conductivity
    property = electrical_conductivity
    block = 'stainless_steel'
  []
  [electrical_conductivity_yttria]
    type = ADMaterialRealAux
    variable = electrical_conductivity
    property = electrical_conductivity
    block = 'powder_compact'
  []
[]

[BCs]
  [external_surface_stainless]
    type = CoupledVarNeumannBC
    boundary = 'outer_radiative_stainless_steel'
    variable = temperature_stainless_steel
    v = heat_transfer_radiation
  []
  [external_surface_graphite]
    type = CoupledVarNeumannBC
    boundary = 'outer_radiative_spacers outer_die_wall radiative_upper_plunger radiative_lower_plunger'
    variable = temperature
    v = heat_transfer_radiation
  []
  [water_channel]
    type = CoupledConvectiveHeatFluxBC
    boundary = water_channel
    variable = temperature_stainless_steel
    T_infinity = 300 #units of K
    htc = 4725
  []
  [temperature_ram_extremes]
    type = ADDirichletBC
    variable = temperature_stainless_steel
    boundary = 'top_upper_steel_ram bottom_lower_steel_ram'
    value = 293
  []
  [electric_potential_top]
    type = ADFunctionNeumannBC
    variable = potential_stainless_steel
    boundary = top_upper_steel_ram
    ## This will need to be updated to match the Dr. Sinter geometry
    function = 'if(t < 31, (980 / (pi * 0.00155))*((0.141301/4.3625)*t), if(t > 600, 0, 980 / (pi * 0.00155)))' # RMS Current / Cross-sectional Area. Ramping for t < 31s approximately reflects Cincotti et al (DOI: 10.1002/aic.11102) Figure 21(b)
  []
  [electric_potential_bottom]
    type = ADDirichletBC
    variable = potential_stainless_steel
    boundary = bottom_lower_steel_ram
    value = 0
  []
[]

[Functions]
  [graphite_thermal_conductivity_fcn]
    type = ParsedFunction
    expression = '-4.418e-12*t^4+2.904e-8*t^3-4.688e-5*t^2-0.0316*t+119.659'
  []
  # [yttria_thermal_conductivity_fcn] #from the multiapp
  #   type = ParsedFunction
  #   expression = '3214.46/(t-147.73)'
  # []
  [harmonic_mean_thermal_conductivity]
    type = ParsedFunction
    expression = '2*(-4.418e-12*t^4+2.904e-8*t^3-4.688e-5*t^2-0.0316*t+119.659)*(3214.46/(t-147.73))/((-4.418e-12*t^4+2.904e-8*t^3-4.688e-5*t^2-0.0316*t+119.659)+(3214.46/(t-147.73)))'
    # symbol_names = 'k_graphite k_yttria'
    # symbol_values = 'graphite_thermal_conductivity_fcn yttria_thermal_conductivity_fcn'
  []

  [graphite_electrical_conductivity_fcn]
    type = ParsedFunction
    expression = '1.0/(-2.705e-15*t^3+1.263e-11*t^2-1.836e-8*t+1.813e-5)'
  []
  # [electrical_conductivity_fcn]
  #   type = ParsedFunction
  #   # symbol_names = porosity
  #   # symbol_values = initial_porosity
  #   expression = '(1-0.62)*2.0025e4*exp(-1.61/8.617343e-5/t)'
  # []
  [harmonic_mean_electrical_conductivity]
    type = ParsedFunction
    expression = '2*(1.0/(-2.705e-15*t^3+1.263e-11*t^2-1.836e-8*t+1.813e-5))*((1)*2.0025e4*exp(-1.61/8.617343e-5/t))/((1.0/(-2.705e-15*t^3+1.263e-11*t^2-1.836e-8*t+1.813e-5))+((1)*2.0025e4*exp(-1.61/8.617343e-5/t)))'
    # symbol_names = 'k_graphite k_yttria'
    # symbol_values = 'graphite_thermal_conductivity_fcn yttria_thermal_conductivity_fcn'
  []
  [mechanical_pressure_func]
    type = ParsedFunction
    symbol_names = 'radius coolant_radius force'
    symbol_values = '0.04 7.071e-3 6.5' # 'm kN'
    expression = 'force * 1e3 / (pi * (radius^2 - coolant_radius^2))' # (N / m^2)
  []
[]

[InterfaceKernels]
  [electric_contact_conductance_ssg]
    type = ElectrostaticContactCondition
    variable = potential_stainless_steel
    neighbor_var = electric_potential
    boundary = 'upper_ram_spacer_interface lower_ram_spacer_interface'
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
  []
  [thermal_contact_conductance_calculated_ssg]
    type = ThermalContactCondition
    variable = temperature_stainless_steel
    neighbor_var = temperature
    primary_potential = potential_stainless_steel
    secondary_potential = electric_potential
    primary_thermal_conductivity = thermal_conductivity
    secondary_thermal_conductivity = thermal_conductivity
    primary_electrical_conductivity = electrical_conductivity
    secondary_electrical_conductivity = electrical_conductivity
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
    boundary = 'upper_ram_spacer_interface lower_ram_spacer_interface'
  []
  [electric_contact_conductance_gss]
    type = ElectrostaticContactCondition
    variable = electric_potential
    neighbor_var = potential_stainless_steel
    boundary = 'upper_ram_spacer_interface lower_ram_spacer_interface'
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
  []
  [thermal_contact_conductance_calculated_gss]
    type = ThermalContactCondition
    variable = temperature
    neighbor_var = temperature_stainless_steel
    primary_potential = electric_potential
    secondary_potential = potential_stainless_steel
    primary_thermal_conductivity = thermal_conductivity
    secondary_thermal_conductivity = thermal_conductivity
    primary_electrical_conductivity = electrical_conductivity
    secondary_electrical_conductivity = electrical_conductivity
    mean_hardness = graphite_stainless_mean_hardness
    mechanical_pressure = mechanical_pressure_func
    boundary = 'upper_ram_spacer_interface lower_ram_spacer_interface'
  []
[]

[ThermalContact]
  [upper_plunger_powder_electric]
    type = GapHeatTransfer
    primary = bottom_upper_plunger
    secondary = top_powder_compact
    variable = electric_potential
    quadrature = true
    emissivity_primary = 0.0 #not applicable for electric potential
    emissivity_secondary = 0.0 #not applicable for electric potential
    gap_geometry_type = PLATE
    # gap_conductivity = 1e-23 #rough harmonic mean
    gap_conductivity_function = 'harmonic_mean_electrical_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [powder_bottom_plunger_electric]
    type = GapHeatTransfer
    primary = top_lower_plunger
    secondary = bottom_powder_compact
    variable = electric_potential
    quadrature = true
    emissivity_primary = 0.0 #not applicable for electric potential
    emissivity_secondary = 0.0 #not applicable for electric potential
    gap_geometry_type = PLATE
    # gap_conductivity = 1e-23 #rough harmonic mean
    gap_conductivity_function = 'harmonic_mean_electrical_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [powder_die_electric]
    type = GapHeatTransfer
    primary = inner_die_wall
    secondary = outer_powder_compact
    variable = electric_potential
    quadrature = true
    emissivity_primary = 0.0 #0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.0 #0.85
    # gap_geometry_type = PLATE # Not for vertical surfaces
    # gap_conductivity = 1e-23 #rough harmonic mean
    gap_conductivity_function = 'harmonic_mean_electrical_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [combined_plungers_die_electric]
    type = GapHeatTransfer
    primary = inner_die_wall
    secondary = combined_outer_plungers
    variable = electric_potential
    quadrature = true
    emissivity_primary = 0.0 #not applicable for electric potential
    emissivity_secondary = 0.0 #not applicable for electric potential
    # gap_geometry_type = PLATE # Not for vertical surfaces
    # gap_conductivity = 8.5e4 #graphite, at 500K
    gap_conductivity_function = 'graphite_electrical_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []

##Thermal Contact between gapped graphite die components
  [upper_plunger_spacer_gap_thermal]
    type = GapHeatTransfer
    primary = spacer_facing_upper_plunger
    secondary = plunger_facing_upper_spacer
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    gap_geometry_type = PLATE
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [upper_diewall_spacer_gap_thermal]
    type = GapHeatTransfer
    primary = top_die_wall
    secondary = bottom_upper_spacer
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    gap_geometry_type = PLATE
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [lower_plunger_spacer_gap_thermal]
    type = GapHeatTransfer
    primary = spacer_facing_lower_plunger
    secondary = plunger_facing_lower_spacer
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    gap_geometry_type = PLATE
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [lower_diewall_spacer_gap_thermal]
    type = GapHeatTransfer
    primary = bottom_die_wall
    secondary = top_lower_spacer
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    gap_geometry_type = PLATE
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [upper_plunger_diewall_gap_thermal]
    type = GapHeatTransfer
    primary = inner_die_wall ### paired temperature doesn't show on inner die wall, but temperature profile looks reasonable
    secondary = die_wall_facing_upper_plunger
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    # gap_geometry_type = PLATE # Not for vertical surfaces
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [lower_plunger_diewall_gap_thermal]
    type = GapHeatTransfer
    primary = inner_die_wall ### paired temperature doesn't show on inner die wall, but temperature profile looks reasonable
    secondary = die_wall_facing_lower_plunger
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    # gap_geometry_type = PLATE # Not for vertical surfaces
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []

## Thermal Contact between touching components of powder and die
  [upper_plunger_powder_thermal]
    type = GapHeatTransfer
    primary = bottom_upper_plunger
    secondary = top_powder_compact
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.3 #estimated from McMahon and Wilder, High Temperature Spectral Emissivity of Yttrium, Samarium, Gadolinium, Ebrium and Lutetium Oxides (1963) Atomic Energy Commission, IS-578, Figure 12
    gap_geometry_type = PLATE
    gap_conductivity_function = 'harmonic_mean_thermal_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [powder_bottom_plunger_thermal]
    type = GapHeatTransfer
    primary = top_lower_plunger
    secondary = bottom_powder_compact #expect more heat transfer from the die to the powder
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.3 #estimated from McMahon and Wilder, High Temperature Spectral Emissivity of Yttrium, Samarium, Gadolinium, Ebrium and Lutetium Oxides (1963) Atomic Energy Commission, IS-578, Figure 12
    gap_geometry_type = PLATE
    gap_conductivity_function = 'harmonic_mean_thermal_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [powder_die_thermal]
    type = GapHeatTransfer
    primary = inner_die_wall
    secondary = outer_powder_compact
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.3 #estimated from McMahon and Wilder, High Temperature Spectral Emissivity of Yttrium, Samarium, Gadolinium, Ebrium and Lutetium Oxides (1963) Atomic Energy Commission, IS-578, Figure 12
    # gap_geometry_type = PLATE # Not for vertical surfaces
    gap_conductivity_function = 'harmonic_mean_thermal_conductivity'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
  [combined_plungers_die_thermal]
    type = GapHeatTransfer
    primary = inner_die_wall
    secondary = combined_outer_plungers
    variable = temperature
    quadrature = true
    emissivity_primary = 0.85 #cincotti 2007, table 2
    emissivity_secondary = 0.85
    # gap_geometry_type = PLATE # Not for vertical surfaces
    gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
    gap_conductivity_function_variable = temperature
    normal_smoothing_distance = 0.1
  []
[]

[Materials]
  ## graphite blocks
  [graphite_thermal]
    type = ADGraphiteThermal
    temperature = temperature
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [graphite_density]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = 1.750e3 #in kg/m^3 from Cincotti et al 2007, Table 2, doi:10.1002/aic
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [graphite_electrical_conductivity]
    type = ADGraphiteElectricalConductivity
    temperature = temperature
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []

  #stainless_steel
  [stainless_steel_thermal]
    type = ADStainlessSteelThermal
    temperature = temperature_stainless_steel
    block = stainless_steel
  []
  [stainless_steel_density]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = 8000
    block = stainless_steel
  []
  [stainless_steel_electrical_conductivity]
    type = ADStainlessSteelElectricalConductivity
    temperature = temperature_stainless_steel
    block = stainless_steel
  []

  # harmonic mean of graphite and stainless steel hardness
  [mean_hardness_graphite_stainless_steel]
    type = GraphiteStainlessMeanHardness
    block = 'graphite_spacers upper_plunger lower_plunger die_wall stainless_steel'
  []

  ## yttria powder compact
  [yttria_thermal_conductivity]
    type = ADParsedMaterial
    coupled_variables = 'temperature'
    expression = '3214.46 / (temperature - 147.73)' #in W/(m-K) #Given from Larry's curve fitting, data from Klein and Croft, JAP, v. 38, p. 1603 and UC report "For Computer Heat Conduction Calculations - A compilation of thermal properties data" by A.L. Edwards, UCRL-50589 (1969)
    # coupled_variables = 'thermal_conductivity_aeh'
    # expression = 'thermal_conductivity_aeh' #in W/(m-K) directly, for now
    property_name = 'thermal_conductivity'
    output_properties = thermal_conductivity
    outputs = 'exodus'
    block = powder_compact
  []
  [yttria_specific_heat_capacity]
    type = ADParsedMaterial
    property_name = heat_capacity
    coupled_variables = 'yttria_heat_capacity_volume_avg'
    expression = 'yttria_heat_capacity_volume_avg' #in J/(K-kg)
    # output_properties = yttria_specific_heat_capacity
    # outputs = 'exodus'
    block = powder_compact
  []
  [yttria_density]
    type = ADParsedMaterial
    property_name = 'density'
    coupled_variables = 'yttria_density_volume_avg'
    expression = 'yttria_density_volume_avg'
    # output_properties = yttria_density
    # outputs = 'exodus'
    block = powder_compact
  []
  [electrical_conductivity]
    type = ADParsedMaterial
    #   coupled_variables = 'yttria_sigma_aeh'
    #   expression = 'yttria_sigma_aeh*1.602e8' #converts to units of J/(V^2-m-s)
    property_name = 'electrical_conductivity'
    output_properties = electrical_conductivity
    outputs = 'exodus'
    block = powder_compact
    # type = ADDerivativeParsedMaterial
    # property_name = electrical_conductivity
    coupled_variables = 'temperature'
    constant_names = 'Q_elec  kB            prefactor_solid  initial_porosity'
    constant_expressions = '1.61    8.617343e-5        1.25e-4           0.38'
    expression = '(1-initial_porosity) * prefactor_solid * exp(-Q_elec/kB/temperature) * 1.602e8' # in eV/(nV^2 s nm) per chat with Larry, last term converts to units of J/(V^2-m-s)
  []

  # Material property converter for DiffusionFluxAux object
  [converter]
    type = MaterialADConverter
    ad_props_in = thermal_conductivity
    reg_props_out = nonad_thermal_conductivity
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  automatic_scaling = true
  line_search = 'none'
  # compute_scaling_once = false

  # force running options
  # petsc_options_iname = '-pc_type -snes_linesearch_type -pc_factor_shift_type -pc_factor_shift_amount'
  # petsc_options_value = 'lu       basic                 NONZERO               1e-15'

  # #mechanical contact options
  # petsc_options = '-snes_ksp_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = ' lu       superlu_dist'
  petsc_options = '-snes_converged_reason -ksp_converged_reason'

  nl_forced_its = 1
  nl_rel_tol = 2e-5 #1e-6 #2e-5 for with mechanics #was 1e-10, for temperature only
  nl_abs_tol = 2e-12 #was 1e-12
  nl_max_its = 20
  l_max_its = 50
  dtmin = 1.0e-4

  end_time = 900 #600 #900 #15 minutes, rule of thumb from Dennis is 10 minutes
  [Quadrature]
    order = FIFTH #required for thermal and mechanical node-face contact
    side_order = SEVENTH
  []
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 0.05
    optimal_iterations = 8
    iteration_window = 2
  []
[]

[Postprocessors]
  [temperature_pp]
    type = AverageNodalVariableValue
    variable = temperature
    block = 'graphite_spacers upper_plunger lower_plunger die_wall powder_compact'
  []
  [temperature_stainless_steel_pp]
    type = AverageNodalVariableValue
    variable = temperature_stainless_steel
    block = stainless_steel
  []
  [graphite_thermal_conductivity]
    type = ElementAverageValue
    variable = thermal_conductivity
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [graphite_electrical_conductivity]
    type = ElementAverageValue
    variable = electrical_conductivity
    block = 'graphite_spacers upper_plunger lower_plunger die_wall'
  []
  [steel_thermal_conductivity]
    type = ElementAverageValue
    variable = thermal_conductivity
    block = stainless_steel
  []
  [steel_electrical_conductivity]
    type = ElementAverageValue
    variable = electrical_conductivity
    block = stainless_steel
  []
  [yttria_thermal_conductivity]
    type = ElementAverageValue
    variable = thermal_conductivity
    block = powder_compact
  []
  [yttria_electrical_conductivity]
    type = ElementAverageValue
    variable = electrical_conductivity
    block = powder_compact
  []

  [yttria_grad_potential]
    type = ElementalVariableValue
    variable = E_y
    elementid = 117
  []
[]

[MultiApps]
  [micro]
    type = TransientMultiApp
    # type = CentroidMultiApp # lauches one in the middle of each element so don't need to give positions
    #can specify the number of procs
    max_procs_per_app = 1 #paolo recommends starting here
    app_type = MalamuteApp
    positions = '0.00125 0.034 0' #roughly the center of element 117 in this mesh
    #positions = '0.0074 0.0058 0' #roughly the center of element 368 in this mesh
    input_files = micro_yttria_thermoelectric_oneway_controls.i
    sub_cycling = true
    execute_on = TIMESTEP_BEGIN #the default
  []
[]

[Transfers]
  [temperature_to_sub]
    type = MultiAppVariableValueSampleTransfer
    to_multi_app = micro
    source_variable = temperature
    variable = T
  []
  [micro_field_pp_to_sub]
    type = MultiAppVariableValueSamplePostprocessorTransfer
    to_multi_app = micro
    source_variable = E_y
    postprocessor = Ey_in
  []
[]

[Outputs]
  csv = true
  exodus = true
  perf_graph = true
  # [ckpt]
  #   type =Checkpoint
  #   time_step_interval = 1
  #   num_files = 2
  # []
[]

Description
Example Input File Syntax
Input Parameters
Input Files
References